Spring 2023 CS4641/CS7641 A Homework 3¶

Instructor: Dr. Mahdi Roozbahani¶

Deadline: Friday, April 7th, 11:59 pm EST¶

  • No unapproved extension of the deadline is allowed. Submission past our 48-hour penalized acceptance period will lead to 0 credit.

  • Discussion is encouraged on Ed as part of the Q/A. However, all assignments should be done individually.

  • Plagiarism is a serious offense. You are responsible for completing your own work. You are not allowed to copy and paste, or paraphrase, or submit materials created or published by others, as if you created the materials. All materials submitted must be your own.</font>
  • All incidents of suspected dishonesty, plagiarism, or violations of the Georgia Tech Honor Code will be subject to the institute’s Academic Integrity procedures. If we observe any (even small) similarities/plagiarisms detected by Gradescope or our TAs, WE WILL DIRECTLY REPORT ALL CASES TO OSI, which may, unfortunately, lead to a very harsh outcome. Consequences can be severe, e.g., academic probation or dismissal, grade penalties, a 0 grade for assignments concerned, and prohibition from withdrawing from the class. </font>

Instructions for the assignment¶

  • This assignment consists of both programming and theory questions.

  • Unless a theory question explicitly states that no work is required to be shown, you must provide an explanation, justification, or calculation for your answer.

  • To switch between cell for code and for markdown, see the menu -> Cell -> Cell Type

  • You can directly type Latex equations into markdown cells.

  • If a question requires a picture, you could use this syntax <img src="" style="width: 300px;"/> to include them within your ipython notebook.

  • Your write up must be submitted in PDF form. You may use either Latex, markdown, or any word processing software. We will **NOT** accept handwritten work. Make sure that your work is formatted correctly, for example submit ∑i=0xi∑i=0xi instead of \text{sum_{i=0} x_i}

  • When submitting the non-programming part of your assignment, you must correctly map pages of your PDF to each question/subquestion to reflect where they appear. **Improperly mapped questions may not be graded correctly and/or will result in point deductions for the error.**
  • All assignments should be done individually, and each student must write up and submit their own answers.
  • Graduate Students: You are required to complete any sections marked as Bonus for Undergrads

Using the autograder¶

  • Grads will find three assignments on Gradescope that correspond to HW3: "Assignment 3 Programming", "Assignment 3 - Non-programming" and "Assignment 3 Programming - Bonus for all". Undergrads will find an additional assignment called "Assignment 3 Programming - Bonus for Undergrads".
  • You will submit your code for the autograder in the Assignment 3 Programming sections. Please refer to the Deliverables and Point Distribution section for what parts are considered required, bonus for undergrads, and bonus for all.

  • We provided you different .py files and we added libraries in those files please DO NOT remove those lines and add your code after those lines. Note that these are the only allowed libraries that you can use for the homework.

  • You are allowed to make as many submissions until the deadline as you like. Additionally, note that the autograder tests each function separately, therefore it can serve as a useful tool to help you debug your code if you are not sure of what part of your implementation might have an issue.

  • For the "Assignment 3 - Non-programming" part, you will need to submit to Gradescope a PDF copy of your Jupyter Notebook with the cells ran. See this EdStem Post for multiple ways on to convert your .ipynb into a .pdf file. Please refer to the Deliverables and Point Distribution section for an outline of the non-programming questions.

  • When submitting to Gradescope, please make sure to mark the page(s) corresponding to each problem/sub-problem. The pages in the PDF should be of size 8.5" x 11", otherwise there may be a deduction in points for extra long sheets.

Using the local tests ¶

  • For some of the programming questions we have included a local test using a small toy dataset to aid in debugging. The local test sample data and outputs are stored in .py files in the local_tests_folder. The actual local tests are stored in localtests.py.
  • There are no points associated with passing or failing the local tests, you must still pass the autograder to get points.
  • It is possible to fail the local test and pass the autograder since the autograder has a certain allowed error tolerance while the local test allowed error may be smaller. Likewise, passing the local tests does not guarantee passing the autograder.
  • You do not need to pass both local and autograder tests to get points, passing the Gradescope autograder is sufficient for credit.
  • It might be helpful to comment out the tests for functions that have not been completed yet.
  • It is recommended to test the functions as it gets completed instead of completing the whole class and then testing. This may help in isolating errors. Do not solely rely on the local tests, continue to test on the autograder regularly as well.

Deliverables and Points Distribution¶

Q1: Image Compression [30pts]¶

Deliverables: imgcompression.py and printed results¶

  • 1.1 Image Compression [20 pts] - programming

    • svd [4pts]

    • compress [4pts]

    • rebuild_svd [4pts]

    • compression_ratio [4pts]

    • recovered_variance_proportion [4pts]

  • 1.2 Black and White [5 pts] non-programming

  • 1.3 Color Image [5 pts] non-programming

Q2: Understanding PCA [20pts]¶

Deliverables: pca.py and written portion¶

  • 2.1 PCA Implementation [10 pts] - programming

    • fit [5pts]

    • transform [2pts]

    • transform_rv [3pts]

  • 2.2 Visualize [5 pts] programming and non-programming

  • 2.3 PCA Reduced Facemask Dataset Analysis [5 pts] non-programming

  • 2.4 PCA Exploration [0 pts]

Q3: Regression and Regularization [80pts: 50pts + 20pts Bonus for Undergrads + 10pts Bonus for All]¶

Deliverables: regression.py and Written portion¶

  • 3.1 Regression and Regularization Implementations [50pts: 30pts + 20pts Bonus for Undergrad] - programming

    • RMSE [5pts]

    • Construct Poly Features 1D [2pts]

    • Construct Poly Features 2D [3pts]

    • Prediction [5pts]

    • Linear Fit Closed Form [5pts]

    • Ridge Fit Closed Form [5pts]

    • Cross Validation [5pts]

    • Linear Gradient Descent [5pts] Bonus for Undergrad

    • Linear Stochastic Gradient Descent [5pts] Bonus for Undergrad

    • Ridge Gradient Descent [5pts] Bonus for Undergrad

    • Ridge Stochastic Gradient Descent [5pts] Bonus for Undergrad

  • 3.2 About RMSE [3 pts] non-programming

  • 3.3 Testing: General Functions and Linear Regression [5 pts] non-programming

  • 3.4 Testing: Ridge Regression [5 pts] non-programming

  • 3.5 Cross Validation [7 pts] non-programming

  • 3.6 Noisy Input Samples in Linear Regression [10 pts] non-programming BONUS FOR ALL

Q4: Naive Bayes and Logistic Regression [35pts]¶

Deliverables: logistic_regression.py and Written portion¶

  • 4.1 Llama Breed Problem using Naive Bayes [5 pts] non-programming

  • 4.2 News Data Sentiment Classification Using Logistic Regression [30 pts] - programming

    • sigmoid [2 pts]

    • bias_augment [3 pts]

    • predict_probs [5 pts]

    • predict_labels [2 pts]

    • loss [3 pts]

    • gradient [3 pts]

    • accuracy [2 pts]

    • evaluate [5 pts]

    • fit [5 pts]

Q5: Noise in PCA and Linear Regression [15pts]¶

Deliverables: Written portion¶

  • 5.1 Slope Functions [5 pts] non-programming

  • 5.2 Error in Y and Error in X and Y [5 pts] non-programming

  • 5.3 Analysis [5 pts] non-programming

Q6: Feature Reduction.py [25pts Bonus for All]¶

Deliverables: feature_reduction.py and Written portion¶

  • 6.1 Feature Reduction [18 pts] - programming

    • forward_selection [9pts]

    • backward_elimination [9pts]

  • 6.2 Feature Selection - Discussion [7 pts] non-programming

Q7: Movie Recommendation with SVD [10pts Bonus for All]¶

Deliverables: svd_recommender.py and Written portion¶

  • 7.1 SVD Recommender

    • recommender_svd [5pts]

    • predict [5pts]

  • 7.2 Visualize Movie Vectors [0pts]

0 Set up¶

This notebook is tested under python 3. ., and the corresponding packages can be downloaded from miniconda. You may also want to get yourself familiar with several packages:

  • jupyter notebook
  • numpy
  • matplotlib
  • sklearn

There is also a VS Code and Anaconda Setup Tutorial on Ed under the "Links" category

Please implement the functions that have raise NotImplementedError, and after you finish the coding, please delete or comment out raise NotImplementedError.

Library imports¶

In [1]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
# This is cell which sets up some of the modules you might need 
# Please do not change the cell or import any additional packages. 

import numpy as np
import pandas as pd
import matplotlib
from matplotlib import pyplot as plt
from sklearn.feature_extraction import text
from sklearn.datasets import load_diabetes, load_breast_cancer, load_iris
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, accuracy_score
import warnings
import sys

print('Version information')

print('python: {}'.format(sys.version))
print('matplotlib: {}'.format(matplotlib.__version__))
print('numpy: {}'.format(np.__version__))

warnings.filterwarnings('ignore')

%matplotlib inline
%load_ext autoreload
%autoreload 2

STUDENT_VERSION = 1
EO_TEXT, EO_FONT, EO_COLOR = 'TA VERSION', 'Arial Black', 'gray', 
EO_ALPHA, EO_SIZE, EO_ROT = 0.7, 90, 40
Version information
python: 3.11.2 | packaged by Anaconda, Inc. | (main, Mar 27 2023, 23:35:04) [MSC v.1916 64 bit (AMD64)]
matplotlib: 3.7.1
numpy: 1.23.5

Q1: Image Compression [30 pts] **[P]** | **[W]**¶

Load images data and plot¶

In [2]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
# load Image
image = plt.imread("./data/hw3_image_compression.jpeg")/255
# plot image
fig = plt.figure(figsize=(10,10))
if not STUDENT_VERSION:
    fig.text(0.5, 0.5, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE, color=EO_COLOR, alpha=EO_ALPHA, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
plt.imshow(image)
Out[2]:
<matplotlib.image.AxesImage at 0x2ab8fb0b790>
In [3]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
def rgb2gray(rgb):   
    return np.dot(rgb[...,:3], [0.299, 0.587, 0.114])

fig = plt.figure(figsize=(10, 10))

if not STUDENT_VERSION:
    fig.text(0.5, 0.5, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE, color=EO_COLOR, alpha=EO_ALPHA, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
# plot several images
plt.imshow(rgb2gray(image), cmap=plt.cm.bone)
Out[3]:
<matplotlib.image.AxesImage at 0x2ab8fba4e10>

1.1 Image compression [20pts] **[P]**¶

SVD is a dimensionality reduction technique that allows us to compress images by throwing away the least important information.

Higher singular values capture greater variance and, thus, capture greater information from the corresponding singular vector. To perform image compression, apply SVD on each matrix and get rid of the small singular values to compress the image. The loss of information through this process is negligible, and the difference between the images can be hardly spotted.

For example, the proportion of variance captured by the first component is σ21∑ni=1σ2iσ12∑i=1nσi2 where σiσi is the ithith singular value.

In the imgcompression.py file, complete the following functions:

  • svd: You may use np.linalg.svd in this function, and although the function defaults this parameter to true, you may explicitly set full_matrices=True using the optional full_matrices parameter. Hint 2 may be useful.
  • compress
  • rebuild_svd
  • compression_ratio: Hint 1 may be useful
  • recovered_variance_proportion: Hint 1 may be useful

HINT 1: http://timbaumann.info/svd-image-compression-demo/ is a useful article on image compression and compression ratio. You may find this article useful for implementing the functions compression_ratio and recovered_variance_proportion

HINT 2: If you have never used np.linalg.svd, it might be helpful to read Numpy's SVD documentation and note the particularities of the VV matrix and that it is returned already transposed.

HINT 3: The shape of SS resulting from SVD may change depending on if N > D, N < D, or N = D. Therefore, when checking the shape of SS, note that min(N,D) means the value should be equal to whichever is lower between N and D.

1.1.1 Local Tests for Image Compression Black and White Case [No Points]¶

You may test your implementation of the functions contained in imgcompression.py in the cell below. Feel free to comment out tests for functions that have not been completed yet. See Using the Local Tests for more details.

In [4]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from utilities.localtests import TestImgCompression

unittest_ic = TestImgCompression()
unittest_ic.test_svd_bw()
unittest_ic.test_compress_bw()
unittest_ic.test_rebuild_svd_bw()
unittest_ic.test_compression_ratio_bw()
unittest_ic.test_recovered_variance_proportion_bw()
UnitTest passed successfully for "SVD calculation - black and white images"!
UnitTest passed successfully for "Image compression - black and white images"!
UnitTest passed successfully for "SVD reconstruction - black and white images"!
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[4], line 11
      9 unittest_ic.test_compress_bw()
     10 unittest_ic.test_rebuild_svd_bw()
---> 11 unittest_ic.test_compression_ratio_bw()
     12 unittest_ic.test_recovered_variance_proportion_bw()

File ~\Desktop\Tech\CS 4641\HW3\hw3_code\utilities\localtests.py:102, in TestImgCompression.test_compression_ratio_bw(self)
     99 test_ic = IC_Test()
    101 # compute compression ratio
--> 102 cr = ic.compression_ratio(test_ic.bw_image, 2)
    104 # check compression ratio
    105 self.assertEqual(
    106     np.allclose(cr, test_ic.cr_g), True, "Compression ratio is incorrect"
    107 )

File ~\Desktop\Tech\CS 4641\HW3\hw3_code\imgcompression.py:111, in ImgCompression.compression_ratio(self, X, k)
    107 N, D = X.shape[1], X.shape[2] if X.ndim == 3 else X.shape
    109 original_size = X.size
--> 111 compressed_size = (N * k + k + k * D) * (3 if X.ndim == 3 else 1)
    113 return compressed_size / original_size

TypeError: unsupported operand type(s) for +: 'int' and 'tuple'

1.1.2 Local Tests for Image Compression Color Case [No Points]¶

You may test your implementation of the functions contained in imgcompression.py in the cell below. Feel free to comment out tests for functions that have not been completed yet. See Using the Local Tests for more details.

In [5]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from utilities.localtests import TestImgCompression

unittest_ic = TestImgCompression()

unittest_ic.test_svd_color()
unittest_ic.test_compress_color()
unittest_ic.test_rebuild_svd_color()
unittest_ic.test_compression_ratio_color()
unittest_ic.test_recovered_variance_proportion_color()
UnitTest passed successfully for "SVD calculation - color images"!
UnitTest passed successfully for "Image compression - color images"!
UnitTest passed successfully for "SVD reconstruction - color images"!
UnitTest passed successfully for "Compression ratio - color images"!
UnitTest passed successfully for "Recovered variance proportion - color images"!

1.2.1 Black and white [5 pts] **[W]**¶

This question will use your implementation of the functions from Q1.1 to generate a set of images compressed to different degrees. You can simply run the below cell without making any changes to it, assuming you have implemented the functions in Q1.1.

Make sure these images are displayed when submitting the PDF version of the Juypter notebook as part of the non-programming submission of this assignment.

In [6]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
from imgcompression import ImgCompression

imcompression = ImgCompression()
bw_image = rgb2gray(image)
U, S, V = imcompression.svd(bw_image)
component_num = [1,2,5,10,20,40,80,160,256]

fig = plt.figure(figsize=(18, 18))

# plot several images
i=0
for k in component_num:
    U_compressed, S_compressed, V_compressed = imcompression.compress(U, S, V, k)
    img_rebuild = imcompression.rebuild_svd(U_compressed, S_compressed, V_compressed)
    c = np.around(imcompression.compression_ratio(bw_image, k), 4)
    r = np.around(imcompression.recovered_variance_proportion(S, k), 3)
    ax = fig.add_subplot(3, 3, i + 1, xticks=[], yticks=[])
    ax.imshow(img_rebuild, cmap=plt.cm.bone)
    ax.set_title(f"{k} Components")
    if not STUDENT_VERSION:
        ax.text(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
            fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA, fontname=EO_FONT,
            ha='center', va='center', rotation=EO_ROT)
    ax.set_xlabel(f"Compression: {c},\nRecovered Variance: {r}")
    i = i+1
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[6], line 18
     16 U_compressed, S_compressed, V_compressed = imcompression.compress(U, S, V, k)
     17 img_rebuild = imcompression.rebuild_svd(U_compressed, S_compressed, V_compressed)
---> 18 c = np.around(imcompression.compression_ratio(bw_image, k), 4)
     19 r = np.around(imcompression.recovered_variance_proportion(S, k), 3)
     20 ax = fig.add_subplot(3, 3, i + 1, xticks=[], yticks=[])

File ~\Desktop\Tech\CS 4641\HW3\hw3_code\imgcompression.py:111, in ImgCompression.compression_ratio(self, X, k)
    107 N, D = X.shape[1], X.shape[2] if X.ndim == 3 else X.shape
    109 original_size = X.size
--> 111 compressed_size = (N * k + k + k * D) * (3 if X.ndim == 3 else 1)
    113 return compressed_size / original_size

TypeError: unsupported operand type(s) for +: 'int' and 'tuple'
<Figure size 1800x1800 with 0 Axes>

1.2.2 Black and White Compression Savings [No Points]¶

This question will use your implementation of the functions from Q1.1 to compare the number of bytes required to represent the SVD decomposition for the original image to the compressed image using different degrees of compression. You can simply run the below cell without making any changes to it, assuming you have implemented the functions in Q1.1.

Running this cell is primarily for your own understanding of the compression process.

In [7]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
from imgcompression import ImgCompression

imcompression = ImgCompression()
bw_image = rgb2gray(image)
U, S, V = imcompression.svd(bw_image)

component_num = [1,2,5,10,20,40,80,160,256]

# Compare memory savings for BW image
for k in component_num:
    og_bytes, comp_bytes, savings = imcompression.memory_savings(bw_image, U, S, V, k)
    comp_ratio = (og_bytes/comp_bytes)
    og_bytes = imcompression.nbytes_to_string(og_bytes)
    comp_bytes = imcompression.nbytes_to_string(comp_bytes)
    savings = imcompression.nbytes_to_string(savings)
    print(f"{k} components: Original Image: {og_bytes} -> Compressed Image: {comp_bytes}, Savings: {savings}, Compression Ratio {comp_ratio:.1f}:1")
1 components: Original Image: 12.207 MB -> Compressed Image: 20.32 KB, Savings: 12.187 MB, Compression Ratio 615.1:1
2 components: Original Image: 12.207 MB -> Compressed Image: 40.641 KB, Savings: 12.167 MB, Compression Ratio 307.6:1
5 components: Original Image: 12.207 MB -> Compressed Image: 101.602 KB, Savings: 12.108 MB, Compression Ratio 123.0:1
10 components: Original Image: 12.207 MB -> Compressed Image: 203.203 KB, Savings: 12.009 MB, Compression Ratio 61.5:1
20 components: Original Image: 12.207 MB -> Compressed Image: 406.406 KB, Savings: 11.81 MB, Compression Ratio 30.8:1
40 components: Original Image: 12.207 MB -> Compressed Image: 812.812 KB, Savings: 11.413 MB, Compression Ratio 15.4:1
80 components: Original Image: 12.207 MB -> Compressed Image: 1.588 MB, Savings: 10.62 MB, Compression Ratio 7.7:1
160 components: Original Image: 12.207 MB -> Compressed Image: 3.175 MB, Savings: 9.032 MB, Compression Ratio 3.8:1
256 components: Original Image: 12.207 MB -> Compressed Image: 5.08 MB, Savings: 7.127 MB, Compression Ratio 2.4:1

1.3.1 Color image [5 pts] **[W]**¶

This section will use your implementation of the functions from Q1.1 to generate a set of images compressed to different degrees. You can simply run the below cell without making any changes to it, assuming you have implemented the functions in Q1.1.

Make sure these images are displayed when submitting the PDF version of the Juypter notebook as part of the non-programming submission of this assignment.

NOTE: You might get warning "Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers)." This warning is acceptable since some of the pixels may go above 1.0 while rebuilding. You should see similar images to original even with such clipping.

HINT 1: Make sure your implementation of recovered_variance_proportion returns an array of 3 floats for a color image.
HINT 2: Try performing SVD on the individual color channels and then stack the individual channel UU, SS, VV matrices.
HINT 3: You may need separate implementations for a color or grayscale image in the same function.

In [8]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
from imgcompression import ImgCompression

imcompression = ImgCompression()
image_rolled = np.moveaxis(image, -1,0)
U, S, V = imcompression.svd(image_rolled)

component_num = [1,2,5,10,20,40,80,160,256]

fig = plt.figure(figsize=(18, 18))

# plot several images
i=0
for k in component_num:
    U_compressed, S_compressed, V_compressed = imcompression.compress(U, S, V, k)
    img_rebuild = np.clip(imcompression.rebuild_svd(U_compressed, S_compressed, V_compressed),0,1)
    img_rebuild = np.moveaxis(img_rebuild, 0,-1)
    c = np.around(imcompression.compression_ratio(image_rolled, k), 4)
    r = np.around(imcompression.recovered_variance_proportion(S, k), 3)
    ax = fig.add_subplot(3, 3, i + 1, xticks=[], yticks=[])
    ax.imshow(img_rebuild)
    ax.set_title(f"{k} Components")
    if not STUDENT_VERSION:
        ax.text(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
            fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA, fontname=EO_FONT,
            ha='center', va='center', rotation=EO_ROT)
    ax.set_xlabel(f"Compression: {np.around(c,4)},\nRecovered Variance:  R: {r[0]}  G: {r[1]}  B: {r[2]}")
    i = i+1

1.3.2 Color Compression Savings [No Points]¶

This question will use your implementation of the functions from Q1.1 to compare the number of bytes required to represent the SVD decomposition for the original image to the compressed image using different degrees of compression. You can simply run the below cell without making any changes to it, assuming you have implemented the functions in Q1.1.

Running this cell is primarily for your own understanding of the compression process.

In [9]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
from imgcompression import ImgCompression

imcompression = ImgCompression()
U, S, V = imcompression.svd(image_rolled)

component_num = [1,2,5,10,20,40,80,160,256]

# Compare the memory savings of the color image
i=0
for k in component_num:
    og_bytes, comp_bytes, savings = imcompression.memory_savings(image_rolled, U, S, V, k)
    comp_ratio = (og_bytes/comp_bytes)
    og_bytes = imcompression.nbytes_to_string(og_bytes)
    comp_bytes = imcompression.nbytes_to_string(comp_bytes)
    savings = imcompression.nbytes_to_string(savings)
    print(f"{k} components: Original Image: {og_bytes} -> Compressed Image: {comp_bytes}, Savings: {savings}, Compression Ratio {comp_ratio:.1f}:1")
1 components: Original Image: 36.621 MB -> Compressed Image: 60.961 KB, Savings: 36.562 MB, Compression Ratio 615.1:1
2 components: Original Image: 36.621 MB -> Compressed Image: 121.922 KB, Savings: 36.502 MB, Compression Ratio 307.6:1
5 components: Original Image: 36.621 MB -> Compressed Image: 304.805 KB, Savings: 36.323 MB, Compression Ratio 123.0:1
10 components: Original Image: 36.621 MB -> Compressed Image: 609.609 KB, Savings: 36.026 MB, Compression Ratio 61.5:1
20 components: Original Image: 36.621 MB -> Compressed Image: 1.191 MB, Savings: 35.43 MB, Compression Ratio 30.8:1
40 components: Original Image: 36.621 MB -> Compressed Image: 2.381 MB, Savings: 34.24 MB, Compression Ratio 15.4:1
80 components: Original Image: 36.621 MB -> Compressed Image: 4.763 MB, Savings: 31.859 MB, Compression Ratio 7.7:1
160 components: Original Image: 36.621 MB -> Compressed Image: 9.525 MB, Savings: 27.096 MB, Compression Ratio 3.8:1
256 components: Original Image: 36.621 MB -> Compressed Image: 15.24 MB, Savings: 21.381 MB, Compression Ratio 2.4:1

Q2: Understanding PCA [20 pts] **[P]** | **[W]**¶

Principal Component Analysis (PCA) is another dimensionality reduction technique that reduces dimensions by eliminating small variance eigenvalues and their vectors. With PCA, we center the data first by subtracting the mean of each feature. Each singular value tells us how much of the variance of a matrix (e.g. image) is captured in each component. In this problem, we will investigate how PCA can be used to improve features for regression and classification tasks and how the data itself affects the behavior of PCA.

Here we will implement PCA using Singular Value Decomposition (SVD). Recall from class that in PCA, we project the original matrix XX into new components, each one corresponding to an eigenvector of the covariance matrix XTXXTX. SVD decomposes XX into three matrices UU, SS and VTVT. This also leads us a decomposition of the covariance matrix as: XTX=(USVT)TUSVT=(VSTUT)USVT=VS2VTXTX=(USVT)TUSVT=(VSTUT)USVT=VS2VT

This means two important things for us:

  • The matrix VTVT, often referred to as the right singular vectors of XX, is equivalent to the eigenvectors of XTXXTX.
  • S2S2 is equivalent to the eigenvalues of XTXXTX.

So the first nn-principal components are obtained by projecting XX by the first nn vectors from VTVT. Similarly, S2S2 gives a measure of the variance retained.

2.1 Implementation [10 pts] **[P]**¶

Implement PCA. In the pca.py file, complete the following functions:

  • fit: You may use np.linalg.svd. Set full_matrices=False. Hint 1 may be useful.
  • transform
  • transform_rv: You may find np.cumsum helpful for this function.

Assume a dataset is composed of N datapoints, each of which has D features with D < N. The dimension of our data would be D. However, it is possible that many of these dimensions contain redundant information. Each feature explains part of the variance in our dataset, and some features may explain more variance than others.

HINT 1: Make sure you remember to first center your data by subtracting the mean of each feature.

2.1.1 Local Tests for PCA [No Points]¶

You may test your implementation of the functions contained in pca.py in the cell below. Feel free to comment out tests for functions that have not been completed yet. See Using the Local Tests for more details.

In [10]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from utilities.localtests import TestPCA

unittest_pca = TestPCA()
unittest_pca.test_pca()
unittest_pca.test_transform()
unittest_pca.test_transform_rv()
UnitTest passed successfully for "PCA fit"!
UnitTest passed successfully for "PCA transform"!
UnitTest passed successfully for "PCA transform with recovered variance"!

2.2 Visualize [5 pts] **[W]**¶

PCA is used to transform multivariate data tables into smaller sets so as to observe the hidden trends and variations in the data. It can also be used as a feature extractor for images. Here you will visualize two datasets using PCA, first is the iris dataset and then a dataset of masked and unmasked images.

In the pca.py, complete the following function:

  • visualize: Use your implementation of PCA and reduce the datasets such that they contain only two features. Using Matplotlib's Pyplot, create 2-D scatter plots of the data points using these features. Make sure to differentiate the data points according to their true labels using color.

The datasets have already been loaded for you in the subsequent cells.

NOTE: Here, we won't be testing for accuracy. Even with correct implementations of PCA, the accuracy can differ from the TA solution. That is fine as long as the visualizations come out similar.

Iris Dataset¶

In [11]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Use PCA for visualization of iris dataset

from pca import PCA

iris_data = load_iris(return_X_y=True)

X = iris_data[0]
y = iris_data[1]

fig = plt.figure()
plt.title('Iris Dataset with Dimensionality Reduction')
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
PCA().visualize(X,y,fig)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[11], line 18
     16 plt.xlabel("Feature 1")
     17 plt.ylabel("Feature 2")
---> 18 PCA().visualize(X,y,fig)

File ~\Desktop\Tech\CS 4641\HW3\hw3_code\pca.py:120, in PCA.visualize(self, X, y, fig)
    104 def visualize(self, X: np.ndarray, y: np.ndarray, fig=None) -> None:  # 5 pts
    105     """
    106     Use your PCA implementation to reduce the dataset to only 2 features. You'll need to run PCA on the dataset and then transform it so that the new dataset only has 2 features.
    107     Create a scatter plot of the reduced data set and differentiate points that have different true labels using color.
   (...)
    117     Return: None
    118     """
--> 120     X_new = self.transform(X, K = 2)
    122     if fig is None:
    123         fig = plt.figure()

File ~\Desktop\Tech\CS 4641\HW3\hw3_code\pca.py:55, in PCA.transform(self, data, K)
     39 """
     40 Transform data to reduce the number of features such that final data (X_new) has K features (columns)
     41 Utilize self.U, self.S and self.V that were set in fit() method.
   (...)
     50 Hint: Make sure you remember to first center your data by subtracting the mean of each feature.
     51 """
     53 X_centered = data - np.mean(data, axis = 0)
---> 55 X_new = np.matmul(X_centered, self.V.T[:, :K])
     56 return X_new

AttributeError: 'NoneType' object has no attribute 'T'

2.3 PCA Reduced Facemask Dataset Analysis [5 pts] **[W]**¶

Facemask Dataset¶

The masked and unmasked dataset is made up of grayscale images of human faces facing forward. Half of these images are faces that are completely unmasked, and the remaining images show half of the face covered with an artificially generated face mask. The images have already been preprocessed, they are also reduced to a small size of 64x64 pixels and then reshaped into a feature vector of 4096 pixels. Below is a sample of some of the images in the dataset.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

X = np.load('./data/smallflat_64.npy')
y = np.load('./data/masked_labels.npy').astype('int')
i = 0
fig = plt.figure(figsize=(18, 18))
for idx in [0,1,2,150,151,152]:
    ax = fig.add_subplot(6, 6, i + 1, xticks=[], yticks=[])
    ax.imshow(X[idx].reshape(64, 64), cmap = 'gray')
    m_status = 'Unmasked' if idx < 150 else 'Masked'
    ax.set_title(f"{m_status} Image at i = {idx}")
    i += 1
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
# Use PCA for visualization of masked and unmasked images

X = np.load('./data/smallflat_64.npy')
y = np.load('./data/masked_labels.npy')

fig = plt.figure()
plt.title('Facemask Dataset Visualization with Dimensionality Reduction')
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
PCA().visualize(X,y, fig)
print('*In this plot, the 0 points are unmasked images and the 1 points are masked images.')

What do you think of this 2 dimensional plot, knowing that the original dataset was originally a set of flattened image vectors that had 4096 pixels/features?.

  1. Look at the 2-dimensional plot above. If the facemask dataset that has been reduced to 2 features was fed into a classifier, do you think the classifier would produce high accuracy or low accuracy in comparison to the original dataset which had 4096 pixels/features? Why? You can refer to the 2D visualization made above (One or two sentences will suffice for this question) (3 pts)

    Answer ...

    The classifier will end up with a lower accuracy because it cannot specify the main patterns and variations in the dataset due to a removal of features

  1. Assuming an equal rate of accuracy, what do you think is the main advantage in feeding a classifier a dataset with 2 features vs a dataset with 4096 features? (One sentence will suffice for this question.) (2 pts)

    Answer ...

    Computational Performance would be significantly higher in 2 features rather than 4096 features.

2.4 PCA Exploration [No Points]¶

Note The accuracy can differ from the TA solution and this section is not graded.

Emotion Dataset [No Points]¶

Now you will use PCA on an actual real-world dataset. We will use your implementation of PCA function to reduce the dataset with 99% retained variance and use it to obtain the reduced features. On the reduced dataset, we will use logistic and linear regression to compare results between PCA and non-PCA datasets. Run the following cells to see how PCA works on regression and classification tasks.

The first dataset we will use is an emotion dataset made up of grayscale images of human faces faces that are visibly happy and visibly sad. Note how Accuracy increases after reducing the number of features used.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

X = np.load('./data/emotion_features.npy')
y = np.load('./data/emotion_labels.npy').astype('int')
i = 0
fig = plt.figure(figsize=(18, 18))
for idx in [0,1,2,150,151,152]:
    ax = fig.add_subplot(6, 6, i + 1, xticks=[], yticks=[])
    ax.imshow(X[idx].reshape(64, 64), cmap = 'gray')
    m_status = 'Sad' if idx < 150 else 'Happy'
    ax.set_title(f"{m_status} Image at i = {idx}")
    i += 1
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

X = np.load('./data/emotion_features.npy')
y = np.load('./data/emotion_labels.npy').astype('int')

print("Not Graded - Data shape before PCA ",X.shape)

pca = PCA()
pca.fit(X)

X_pca = pca.transform_rv(X, retained_variance = 0.99)

print("Not Graded - Data shape with PCA ",X_pca.shape)
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
# Train, test splits
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.3, 
                                                    stratify=y, 
                                                    random_state=42)

# Use logistic regression to predict classes for test set
clf = LogisticRegression()
clf.fit(X_train, y_train)
preds = clf.predict_proba(X_test)
print('Not Graded - Accuracy before PCA: {:.5f}'.format(accuracy_score(y_test, 
                                                preds.argmax(axis=1))))
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
# Train, test splits
X_train, X_test, y_train, y_test = train_test_split(X_pca, y, test_size=.3, 
                                                    stratify=y, 
                                                    random_state=42)

# Use logistic regression to predict classes for test set
clf = LogisticRegression()
clf.fit(X_train, y_train)
preds = clf.predict_proba(X_test)
print('Not Graded - Accuracy after PCA: {:.5f}'.format(accuracy_score(y_test, 
                                                preds.argmax(axis=1))))

Now we will explore sklearn's Diabetes dataset using PCA dimensionality reduction and regression. Notice the RMSE score reduction after we apply PCA.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
from sklearn.linear_model import RidgeCV
def apply_regression(X_train, y_train, X_test):
    ridge =  RidgeCV(alphas=[1e-3, 1e-2, 1e-1, 1])
    clf = ridge.fit(X_train, y_train)
    y_pred = ridge.predict(X_test)
    
    return y_pred
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
#load the dataset 
diabetes = load_diabetes()
X = diabetes.data
y = diabetes.target

print(X.shape, y.shape)

pca = PCA()
pca.fit(X)

X_pca = pca.transform_rv(X, retained_variance = 0.9)
print("Not Graded - data shape with PCA ",X_pca.shape)
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
# Train, test splits
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.3, random_state=42)

#Ridge regression without PCA
y_pred = apply_regression(X_train, y_train, X_test)

# calculate RMSE 
rmse_score = np.sqrt(mean_squared_error(y_pred, y_test))
print('Not Graded - RMSE score using Ridge Regression before PCA: {:.5}'.format(rmse_score))
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
#Ridge regression with PCA
X_train, X_test, y_train, y_test = train_test_split(X_pca, y, test_size=.3, random_state=42)

#use Ridge Regression for getting predicted labels
y_pred = apply_regression(X_train,y_train,X_test)

#calculate RMSE 
rmse_score = np.sqrt(mean_squared_error(y_pred, y_test))
print('Not Graded - RMSE score using Ridge Regression after PCA: {:.5}'.format(rmse_score))

Q3 Polynomial regression and regularization [80pts: 50pts + 20pts Bonus for Undergrads + 10pts Bonus for All] **[P]** | **[W]**¶

3.1 Regression and regularization implementations [50pts: 30 pts + 20 pts bonus for CS 4641] **[P]**¶

We have three methods to fit linear and ridge regression models: 1) closed form solution; 2) gradient descent (GD); 3) stochastic gradient descent (SGD). Some of the functions are bonus, see the below function list on what is required to be implemented for graduate and undergraduate students. We use the term weight in the following code. Weights and parameters (θθ) have the same meaning here. We used parameters (θθ) in the lecture slides.

In the regression.py file, complete the Regression class by implementing the listed functions below. We have provided the Loss function, LL, associated with the GD and SGD function for Linear and Ridge Regression for deriving the gradient update.

  • rmse
  • construct_polynomial_feats
  • predict
  • linear_fit_closed: You should use np.linalg.pinv in this function
  • linear_fit_GD (bonus for undergrad, required for grad): Llinear, GD(θ)=12NN∑i=0[yi−^yi(θ)]2yi=label,^yi(θ)=predictionLlinear, GD(θ)=12N∑i=0N[yi−y^i(θ)]2yi=label,y^i(θ)=prediction

  • linear_fit_SGD (bonus for undergrad, required for grad): Llinear, SGD(θ)=12[yi−^yi(θ)]2yi=label,^yi(θ)=predictionLlinear, SGD(θ)=12[yi−y^i(θ)]2yi=label,y^i(θ)=prediction

  • ridge_fit_closed: You should adjust your I matrix to handle the bias term differently than the rest of the terms

  • ridge_fit_GD (bonus for undergrad, required for grad): Lridge, GD(θ)=Llinear, GD(θ)+cλ2NθTθLridge, GD(θ)=Llinear, GD(θ)+cλ2NθTθ

  • ridge_fit_SGD (bonus for undergrad, required for grad):

Lridge, SGD(θ)=Llinear, SGD(θ)+cλ2NθTθLridge, SGD(θ)=Llinear, SGD(θ)+cλ2NθTθ
  • ridge_cross_validation: Use ridge_fit_closed for this function

IMPORTANT NOTE:

  • Use your RMSE function to calculate actual loss when coding GD and SGD, but use the loss listed above to derive the gradient update.

The points for each function is in the Deliverables and Points Distribution section.

3.1.1 Local Tests for Helper Regression Functions [No Points]¶

You may test your implementation of the functions contained in regression.py in the cell below. Feel free to comment out tests for functions that have not been completed yet. See Using the Local Tests for more details.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from utilities.localtests import TestRegression

unittest_reg = TestRegression()
unittest_reg.test_rmse()
unittest_reg.test_construct_polynomial_feats()
unittest_reg.test_predict()

3.1.2 Local Tests for Linear Regression Functions [No Points]¶

You may test your implementation of the functions contained in regression.py in the cell below. Feel free to comment out tests for functions that have not been completed yet. See Using the Local Tests for more details.

In [12]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from utilities.localtests import TestRegression

unittest_reg = TestRegression()
unittest_reg.test_linear_fit_closed()
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[12], line 8
      5 from utilities.localtests import TestRegression
      7 unittest_reg = TestRegression()
----> 8 unittest_reg.test_linear_fit_closed()

File ~\Desktop\Tech\CS 4641\HW3\hw3_code\utilities\localtests.py:477, in TestRegression.test_linear_fit_closed(self)
    470 test_reg = Regression_Test()
    472 linear_closed_test = np.allclose(
    473     reg.linear_fit_closed(test_reg.x_all_feat, test_reg.y_all),
    474     test_reg.linear_closed,
    475     rtol=1e-04,
    476 )
--> 477 self.assertTrue(linear_closed_test, "Weights are incorrect")
    479 success_msg = "Closed form linear regression"
    480 print_success_message(success_msg)

File ~\miniconda3\envs\ml_hw3\Lib\unittest\case.py:715, in TestCase.assertTrue(self, expr, msg)
    713 if not expr:
    714     msg = self._formatMessage(msg, "%s is not true" % safe_repr(expr))
--> 715     raise self.failureException(msg)

AssertionError: False is not true : Weights are incorrect

3.1.3 Local Tests for Ridge Regression Functions [No Points]¶

You may test your implementation of the functions contained in regression.py in the cell below. Feel free to comment out tests for functions that have not been completed yet. See Using the Local Tests for more details.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from utilities.localtests import TestRegression

unittest_reg = TestRegression()
unittest_reg.test_ridge_fit_closed()
unittest_reg.test_ridge_cross_validation()

3.1.4 Local Tests for Gradient Descent and SGD (Bonus for Undergrad Tests) [No Points]¶

You may test your implementation of the functions contained in regression.py in the cell below. Feel free to comment out tests for functions that have not been completed yet. See Using the Local Tests for more details.

In [13]:
from utilities.localtests import TestRegression

unittest_reg = TestRegression()
unittest_reg.test_linear_fit_GD()
unittest_reg.test_linear_fit_SGD()
unittest_reg.test_ridge_fit_GD()
unittest_reg.test_ridge_fit_SGD()
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[13], line 4
      1 from utilities.localtests import TestRegression
      3 unittest_reg = TestRegression()
----> 4 unittest_reg.test_linear_fit_GD()
      5 unittest_reg.test_linear_fit_SGD()
      6 unittest_reg.test_ridge_fit_GD()

File ~\Desktop\Tech\CS 4641\HW3\hw3_code\utilities\localtests.py:496, in TestRegression.test_linear_fit_GD(self)
    494 lgd_loss_test = np.allclose(linear_GD_loss, test_reg.linear_GD_loss)
    495 self.assertTrue(lgd_test, "Weights are incorrect")
--> 496 self.assertTrue(lgd_loss_test, "Loss is incorrect")
    498 success_msg = "Gradient descent linear regression"
    499 print_success_message(success_msg)

File ~\miniconda3\envs\ml_hw3\Lib\unittest\case.py:715, in TestCase.assertTrue(self, expr, msg)
    713 if not expr:
    714     msg = self._formatMessage(msg, "%s is not true" % safe_repr(expr))
--> 715     raise self.failureException(msg)

AssertionError: False is not true : Loss is incorrect

3.2 About RMSE [3 pts] **[W]**¶

What is a good RMSE value?

If we normalize our labels such that the true labels yy and the model outputs ^yy^ can only be between 0 and 1, what does it mean when the RMSE = 1? Please provide an example with your explanation.

**ANSWER:** ...

Having an RMSE of 1 indicates that the model is performing terribly. The RMSE value should be closer to 0 because this indicates that the model is making less errors while a higher RMSE value indicates a larger error in prediction.

3.3 Testing: General Functions and Linear Regression [5 pts] **[W]**¶

In this section. we will test the performance of the linear regression. As long as your test RMSE score is close to the TA's answer (TA's answer ±0.05±0.05), you can get full points. Let's first construct a dataset for polynomial regression.

In this case, we construct the polynomial features up to degree 5. Each data sample consists of two features [a,b][a,b]. We compute the polynomial features of both aa and bb in order to yield the vectors [1,a,a2,a3,…,adegree][1,a,a2,a3,…,adegree] and [1,b,b2,b3,…,bdegree][1,b,b2,b3,…,bdegree]. We train our model with the cartesian product of these polynomial features. The cartesian product generates a new feature vector consisting of all polynomial combinations of the features with degree less than or equal to the specified degree.

For example, if degree = 2, we will have the polynomial features [1,a,a2][1,a,a2] and [1,b,b2][1,b,b2] for the datapoint [a,b][a,b]. The cartesian product of these two vectors will be [1,a,b,ab,a2,b2][1,a,b,ab,a2,b2]. We do not generate a3a3 and b3b3 since their degree is greater than 2 (specified degree).

In [26]:
from regression import Regression
In [27]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
POLY_DEGREE = 7
N_SAMPLES = 1200

rng = np.random.RandomState(seed=10)

# Simulating a regression dataset with polynomial features.
true_weight = rng.rand(POLY_DEGREE ** 2 + 2, 1)
x_feature1 = np.linspace(-5, 5, N_SAMPLES)
x_feature2 = np.linspace(-3, 3, N_SAMPLES)
x_all = np.stack((x_feature1, x_feature2), axis=1)

reg = Regression()
x_all_feat = reg.construct_polynomial_feats(x_all, POLY_DEGREE)
x_cart_flat = []
for i in range(x_all_feat.shape[0]):
    point = x_all_feat[i]
    x1 = point[:,0]
    x2 = point[:,1]
    x1_end = x1[-1]
    x2_end = x2[-1]
    x1 = x1[:-1]
    x2 = x2[:-1]
    x3 = np.asarray([[m*n for m in x1] for n in x2])

    x3_flat = list(np.reshape(x3, (x3.shape[0] ** 2)))
    x3_flat.append(x1_end)
    x3_flat.append(x2_end)
    x3_flat = np.asarray(x3_flat)
    x_cart_flat.append(x3_flat)
  
x_cart_flat = np.asarray(x_cart_flat)
x_cart_flat = (x_cart_flat - np.mean(x_cart_flat)) / np.std(x_cart_flat)  # Normalize
x_all_feat = np.copy(x_cart_flat)

# We must add noise to data, else the data will look unrealistically perfect.
y_noise = rng.randn(x_all_feat.shape[0], 1)
y_all = np.dot(x_cart_flat, true_weight) + y_noise
print("x_all: ", x_all.shape[0], " (rows/samples) ", x_all.shape[1], " (columns/features)", sep="")
print("y_all: ", y_all.shape[0], " (rows/samples) ", y_all.shape[1], " (columns/features)", sep="")
x_all: 1200 (rows/samples) 2 (columns/features)
y_all: 1200 (rows/samples) 1 (columns/features)
In [28]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

p = np.reshape(np.dot(x_cart_flat, true_weight), (N_SAMPLES,))
ax.scatter(x_all[:,0], x_all[:,1], y_all, label='Datapoints', s=4, alpha=0.2)
ax.plot(x_all[:,0], x_all[:,1], p, label='Line of Best Fit', c="red", linewidth=2)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)

ax.legend()
ax.text2D(0.05, 0.95, "All Simulated Datapoints", transform=ax.transAxes)
plt.show()

In the figure above, the red curve is the true fuction we want to learn, while the blue dots are the noisy data points. The data points are generated by Y=Xθ+ϵY=Xθ+ϵ, where ϵi∼N(0,1)ϵi∼N(0,1) are i.i.d. generated noise.

Now let's split the data into two parts, the training set and testing set. The yellow dots are for training, while the black dots are for testing.

In [29]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
PERCENT_TRAIN = 0.8

all_indices = rng.permutation(N_SAMPLES)  # Random indicies
train_indices = all_indices[:round(N_SAMPLES * PERCENT_TRAIN)]  # 80% Training
test_indices = all_indices[round(N_SAMPLES * PERCENT_TRAIN):]  # 20% Testing

xtrain = x_all[train_indices]
ytrain = y_all[train_indices]
xtest = x_all[test_indices]
ytest = y_all[test_indices]

# -- Plotting Code --
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

ax.scatter(xtrain[:,0], xtrain[:,1], ytrain, label='Training', c='y',s=4)
ax.scatter(xtest[:,0], xtest[:,1], ytest, label='Testing', c='black',s=4)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
    
ax.legend(loc = 'upper right')
plt.show()

Now let us train our model using the training set and see how our model performs on the testing set. Observe the red line, which is our model's learned function.

In [30]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Required for both Grad and Undergrad

weight = reg.linear_fit_closed(x_all_feat[train_indices], y_all[train_indices])
y_test_pred = reg.predict(x_all_feat[test_indices], weight)
test_rmse = reg.rmse(y_test_pred, y_all[test_indices])
print('Linear (closed) RMSE: %.4f' % test_rmse)

# -- Plotting Code --
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

y_pred = reg.predict(x_all_feat, weight)
y_pred = np.reshape(y_pred, (y_pred.size,))
ax.plot(x_all[:,0], x_all[:,1], y_pred, label='Trendline', color='r', lw=2, zorder=5)

ax.scatter(xtrain[:,0], xtrain[:,1], ytrain, label='Training', c='y',s=4)
ax.scatter(xtest[:,0], xtest[:,1], ytest, label='Testing', c='black',s=4)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)

ax.text2D(0.05, 0.95, "Linear (Closed)", transform=ax.transAxes)
ax.legend(loc = 'upper right')
plt.show()
Linear (closed) RMSE: 0.9097

HINT: If your RMSE is off, make sure to follow the instruction given for linear_fit_closed in the list of functions to implement above.

Now let's use our linear gradient descent function with the same setup. Observe that the trendline is now less optimal, and our RMSE decreased. Do not be alarmed.

In [31]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Required for Grad Only
# This cell may take more than 1 minute

weight, _ = reg.linear_fit_GD(x_all_feat[train_indices],
                           y_all[train_indices],
                           epochs=50000,
                           learning_rate=1e-8)
y_test_pred = reg.predict(x_all_feat[test_indices], weight)
test_rmse = reg.rmse(y_test_pred, y_all[test_indices])
print('Linear (GD) RMSE: %.4f' % test_rmse)

# -- Plotting Code --
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

y_pred = reg.predict(x_all_feat, weight)
y_pred = np.reshape(y_pred, (y_pred.size,))
ax.plot(x_all[:,0], x_all[:,1], y_pred, label='Trendline', color='r', lw=2, zorder=5)

ax.scatter(xtrain[:,0], xtrain[:,1], ytrain, label='Training', c='y',s=4)
ax.scatter(xtest[:,0], xtest[:,1], ytest, label='Testing', c='black',s=4)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)

ax.text2D(0.05, 0.95, "Linear (GD)", transform=ax.transAxes)
ax.legend(loc = 'upper right')
plt.show()
Linear (GD) RMSE: 5.1484

We must tune our epochs and learning_rate. As we tune these parameters our trendline will approach the trendline generated by the linear closed form solution. Observe how we slowly tune (increase) the epochs and learning_rate below to create a better model.

Note that the closed form solution will always give the most optimal/overfit results. We cannot outperform the closed form solution with GD. We can only approach closed forms level of optimality/overfitness. We leave the reasoning behind this as an exercise to the reader.

In [32]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Required for Grad Only
# This cell may take more than 1 minute

learning_rates = [1e-8, 1e-6, 1e-4]
weights = np.zeros((3, POLY_DEGREE ** 2 + 2))

for ii in range(len(learning_rates)):
    weights[ii,:] = reg.linear_fit_GD(x_all_feat[train_indices],
                                      y_all[train_indices],
                                      epochs=50000,
                                      learning_rate=learning_rates[ii])[0].ravel()
    y_test_pred = reg.predict(x_all_feat[test_indices],
                              weights[ii, :].reshape((POLY_DEGREE ** 2 + 2, 1)))
    test_rmse = reg.rmse(y_test_pred, y_all[test_indices])
    print('Linear (GD) RMSE: %.4f (learning_rate=%s)' % (test_rmse, learning_rates[ii]))

# -- Plotting Code --
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

colors = ['g', 'orange', 'r']
for ii in range(len(learning_rates)):
    y_pred = reg.predict(x_all_feat, weights[ii])
    y_pred = np.reshape(y_pred, (y_pred.size,))
    ax.plot(x_all[:,0], x_all[:,1], y_pred,
            label='Trendline LR=' + str(learning_rates[ii]),
            color=colors[ii], lw=2, zorder=5)

ax.scatter(xtrain[:,0], xtrain[:,1], ytrain, label='Training', c='y',s=4)
ax.scatter(xtest[:,0], xtest[:,1], ytest, label='Testing', c='black',s=4)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
    
ax.text2D(0.05, 0.95, "Tuning Linear (GD)", transform=ax.transAxes)
ax.legend(loc = 'upper right')
plt.show()
Linear (GD) RMSE: 5.1484 (learning_rate=1e-08)
Linear (GD) RMSE: 3.3447 (learning_rate=1e-06)
Linear (GD) RMSE: 1.1079 (learning_rate=0.0001)

And what if we just use the first 10 data points to train?

In [33]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
rng = np.random.RandomState(seed=5)
y_all_noisy = np.dot(x_cart_flat, np.zeros((POLY_DEGREE ** 2 + 2, 1))) + rng.randn(x_all_feat.shape[0], 1)
sub_train = train_indices[10:20]
In [34]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Required for both Grad and Undergrad

weight = reg.linear_fit_closed(x_all_feat[sub_train], y_all_noisy[sub_train])
y_pred = reg.predict(x_all_feat, weight)
y_test_pred = reg.predict(x_all_feat[test_indices], weight)
test_rmse = reg.rmse(y_test_pred, y_all_noisy[test_indices])
print('Linear (closed) 10 Samples RMSE: %.4f' % test_rmse)

# -- Plotting Code --
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

x1 = x_all[:,0]
x2 = x_all[:,1]
y_pred = np.reshape(y_pred, (N_SAMPLES,))
ax.plot(x1, x2, y_pred, color='b', lw=4)

x3 = x_all[sub_train,0]
x4 = x_all[sub_train,1]
ax.scatter(x3, x4, y_all_noisy[sub_train], s=100, c='r', marker='x')

y_test_pred = reg.predict(x_all_feat[test_indices], weight)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")
ax.set_zlim([None, 8])

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)

ax.text2D(0.05, 0.95, "Linear Regression (Closed)", transform=ax.transAxes)
Linear (closed) 10 Samples RMSE: 806.8299
Out[34]:
Text(0.05, 0.95, 'Linear Regression (Closed)')

Did you see a worse performance? Let's take a closer look at what we have learned.

3.4 Testing: Testing ridge regression [5 pts] **[W]**¶

Now let's try ridge regression. Like before, undergraduate students need to implement the closed form, and graduate students need to implement all three methods. We will call the prediction function from linear regression part. As long as your test RMSE score is close to the TA's answer (TA's answer ±0.05±0.05), you can get full points.

Again, let's see what we have learned. You only need to run the cell corresponding to your specific implementation.

In [35]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Required for both Grad and Undergrad

weight = reg.ridge_fit_closed(x_all_feat[sub_train],
                              y_all_noisy[sub_train],
                              c_lambda=10)
y_pred = reg.predict(x_all_feat, weight)
y_test_pred = reg.predict(x_all_feat[test_indices], weight)
test_rmse = reg.rmse(y_test_pred, y_all_noisy[test_indices])
print('Ridge Regression (closed) RMSE: %.4f' % test_rmse)

# -- Plotting Code --
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

x1 = x_all[:,0]
x2 = x_all[:,1]
y_pred = np.reshape(y_pred, (N_SAMPLES,))
ax.plot(x1, x2, y_pred, color='b', lw=4)

x3 = x_all[sub_train,0]
x4 = x_all[sub_train,1]
ax.scatter(x3, x4, y_all_noisy[sub_train], s=100, c='r', marker='x')

y_test_pred = reg.predict(x_all_feat[test_indices], weight)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
ax.set_zlim([None, 8])
ax.text2D(0.05, 0.95, "Ridge Regression (Closed)", transform=ax.transAxes)
Ridge Regression (closed) RMSE: 1.4765
Out[35]:
Text(0.05, 0.95, 'Ridge Regression (Closed)')

HINT: Make sure to follow the instruction given for ridge_fit_closed in the list of functions to implement above.

In [36]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Required for Grad Only

weight, _ = reg.ridge_fit_GD(x_all_feat[sub_train],
                          y_all_noisy[sub_train],
                          c_lambda=20, learning_rate=1e-5)
y_pred = reg.predict(x_all_feat, weight)
y_test_pred = reg.predict(x_all_feat[test_indices], weight)
test_rmse = reg.rmse(y_test_pred, y_all_noisy[test_indices])
print('Ridge Regression (GD) RMSE: %.4f' % test_rmse)

# -- Plotting Code --
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

x1 = x_all[:,0]
x2 = x_all[:,1]
y_pred = np.reshape(y_pred, (N_SAMPLES,))
ax.plot(x1, x2, y_pred, color='b', lw=4)

x3 = x_all[sub_train,0]
x4 = x_all[sub_train,1]
ax.scatter(x3, x4, y_all_noisy[sub_train], s=100, c='r', marker='x')

y_test_pred = reg.predict(x_all_feat[test_indices], weight)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")
ax.set_zlim([None, 8])

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)

ax.text2D(0.05, 0.95, "Ridge Regression (GD)", transform=ax.transAxes)
Ridge Regression (GD) RMSE: 0.9422
Out[36]:
Text(0.05, 0.95, 'Ridge Regression (GD)')
In [49]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Required for Grad Only

weight, _ = reg.ridge_fit_SGD(x_all_feat[sub_train],
                           y_all_noisy[sub_train],
                           c_lambda=20,
                           learning_rate=1e-5)
y_pred = reg.predict(x_all_feat, weight)
y_test_pred = reg.predict(x_all_feat[test_indices], weight)
test_rmse = reg.rmse(y_test_pred, y_all_noisy[test_indices])
print('Ridge Regression (SGD) RMSE: %.4f' % test_rmse)


# -- Plotting Code --
fig = plt.figure(figsize=(8,5), dpi=120)
ax = fig.add_subplot(111, projection='3d')

x1 = x_all[:,0]
x2 = x_all[:,1]
y_pred = np.reshape(y_pred, (N_SAMPLES,))
ax.plot(x1, x2, y_pred, color='b', lw=4)

x3 = x_all[sub_train,0]
x4 = x_all[sub_train,1]
ax.scatter(x3, x4, y_all_noisy[sub_train], s=100, c='r', marker='x')

y_test_pred = reg.predict(x_all_feat[test_indices], weight)
ax.set_xlabel("feature 1")
ax.set_ylabel("feature 2")
ax.set_zlabel("y")
ax.set_zlim([None, 8])

if not STUDENT_VERSION:
    ax.text2D(0.5, 0.5, EO_TEXT, transform=ax.transAxes,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.4, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
    
ax.text2D(0.05, 0.95, "Ridge Regression (SGD)", transform=ax.transAxes)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[49], line 7
      1 ###############################
      2 ### DO NOT CHANGE THIS CELL ###
      3 ###############################
      4 
      5 # Required for Grad Only
----> 7 weight, _ = reg.ridge_fit_SGD(x_all_feat[sub_train],
      8                            y_all_noisy[sub_train],
      9                            c_lambda=20,
     10                            learning_rate=1e-5)
     11 y_pred = reg.predict(x_all_feat, weight)
     12 y_test_pred = reg.predict(x_all_feat[test_indices], weight)

File ~\Desktop\Tech\CS 4641\HW3\hw3_code\regression.py:310, in Regression.ridge_fit_SGD(self, xtrain, ytrain, c_lambda, epochs, learning_rate)
    308         loss = self.rmse(pred, ytrain[i])
    309         loss_per_step.append(loss)
--> 310         weight = weight * (1 - learning_rate * c_lambda / N) - learning_rate * xtrain[i].T @ (pred - ytrain[i]) / N
    312 return weight, loss_per_step

ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 1 is different from 51)

3.5 Cross validation [7 pts] **[W]**¶

Let's use Cross Validation to search for the best value for c_lambda in ridge regression.

Imagine we have a dataset of 10 points [1,2,3,4,5,6,7,8,9,10] and we want to do 5-fold cross validation.

  • The first iteration we would train with [3,4,5,6,7,8,9,10] and test (validate) with [1,2]
  • The second iteration we would train with [1,2,5,6,7,8,9,10] and test (validate) with [3,4]
  • The third iteration we would train with [1,2,3,4,7,8,9,10] and test (validate) with [5,6]
  • The fourth iteration we would train with [1,2,3,4,5,6,9,10] and test (validate) with [7,8]
  • The fifth iteration we would train with [1,2,3,4,5,6,7,8] and test (validate) with [9,10]

We provided a list of possible values for λλ, and you will use them in cross validation. For cross validation, use 10-fold method and only use it for your training data (you already have the train_indices to get training data). For the training data, split them in 10 folds which means that use 10 percent of training data for test and 90 percent for training. For each λλ, you will have calculated 10 RMSE values. Compute the mean of the 10 RMSE values. Then pick the λλ with the lowest mean RMSE.

HINTS:

  • np.concatenate is your friend
  • Make sure to follow the instruction given for ridge_fit_closed in the list of functions to implement above.
  • To use the 10-fold method, that would include looping over all the data 10 times, where we split a different 10% of the data at every iteration. So the first iteration extracts the first 10% to testing and the remaining 90% for training.The second iteration splits the second 10% of data for testing and the (different) remaining 90% for testing. If we have the array of elements 1 - 10, the second iteration would extract the number "2" because that's in the second 10% of the array.
  • The hyperparameter_search function will handle averaging the errors, so don't average the errors in ridge_cross_validation. We've done this so you can see your error across every fold when using the gradescope tests.
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
lambda_list = [0.0001, 0.001, 0.1, 1, 5, 10, 50, 100, 1000, 10000]
kfold = 10

best_lambda, best_error, error_list = reg.hyperparameter_search(x_all_feat[train_indices], y_all[train_indices], lambda_list, kfold)
for lm, err in zip(lambda_list, error_list):
    print('Lambda: %.4f' % lm, 'RMSE: %.6f'% err)

print('Best Lambda: %.4f' % best_lambda)
weight = reg.ridge_fit_closed(x_all_feat[train_indices], y_all_noisy[train_indices], c_lambda=best_lambda)
y_test_pred = reg.predict(x_all_feat[test_indices], weight)
test_rmse = reg.rmse(y_test_pred, y_all_noisy[test_indices])
print('Best Test RMSE: %.4f' % test_rmse) 

3.6 Noisy Input Samples in Linear Regression [10 pts Bonus for All] **[W]**¶

Consider a linear model of the form: y(xn,θ)=θ0+D∑d=1θdxndy(xn,θ)=θ0+∑d=1Dθdxnd where xn=(xn1,…,xnD)∈RDxn=(xn1,…,xnD)∈RD and weights θ=(θ0,…,θD)∈RDθ=(θ0,…,θD)∈RD. Given the the D-dimension input sample set x={x1,…,xn}x={x1,…,xn} with corresponding target value y={y1,…,yn}y={y1,…,yn}, the sum-of-squares error function is: ED(θ)=12N∑n=1[y(xn,θ)−yn]2ED(θ)=12∑n=1N[y(xn,θ)−yn]2

Now, suppose that Gaussian noise ϵn∈RDϵn∈RD is added independently to each of the input sample xnxn to generate a new sample set x′={x1+ϵ1,…,xn+ϵn}x′={x1+ϵ1,…,xn+ϵn}. Here, ϵniϵni (an entry of ϵnϵn) has zero mean and variance σ2σ2. For each sample xnxn, let x′n=(xn1+ϵn1,…,xnD+ϵnd)xn′=(xn1+ϵn1,…,xnD+ϵnd), where nn and dd is independent across both nn and dd indices.

  1. (3pts) Show that y(x′n,θ)=y(xn,θ)+∑Dd=1θdϵndy(xn′,θ)=y(xn,θ)+∑d=1Dθdϵnd

  2. (7pts) Assume the sum-of-squares error function of the noise sample set x′={x1+ϵ1,…,xn+ϵn}x′={x1+ϵ1,…,xn+ϵn} is ED(θ)′ED(θ)′. Prove the expectation of ED(θ)′ED(θ)′ is equivalent to the sum-of-squares error ED(θ)ED(θ) for noise-free input samples with the addition of a weight-decay regularization term (e.g. ℓ2ℓ2 norm) , in which the bias parameter θ0θ0 is omitted from the regularizer. In other words, show that E[ED(θ)′]=ED(θ)+Regularizer.E[ED(θ)′]=ED(θ)+Regularizer.

N.B. You should be incorporating your solution from the first part of this problem into the given sum of squares equation for the second part.

HINT:

  • During the class, we have discussed how to solve for the weight θθ for ridge regression, the function looks like this: E(θ)=1NN∑i=1[y(xi,θ)−yi]2+λNd∑i=1|θi|2E(θ)=1N∑i=1N[y(xi,θ)−yi]2+λN∑i=1d|θi|2 where the first term is the sum-of-squares error and the second term is the regularization term. N is the number of samples. In this question, we use another form of the ridge regression, which is: E(θ)=12N∑i=1[y(xi,θ)−yi]2+λ2d∑i=1|θi|2E(θ)=12∑i=1N[y(xi,θ)−yi]2+λ2∑i=1d|θi|2
  • For the Gaussian noise ϵnϵn, we have E[ϵn]=0E[ϵn]=0

  • Assume the noise ϵ=(ϵ1,...,ϵn)ϵ=(ϵ1,...,ϵn) are independent to each other, we have E[ϵnϵm]={σ2m=n0m≠nE[ϵnϵm]={σ2m=n0m≠n

  1. Answer: ...
  1. Answer: ...

Q4: Naive Bayes and Logistic Regression [35pts] **[P]** | **[W]**¶

In Bayesian classification, we're interested in finding the probability of a label given some observed feature vector x=[x1,…,xd]x=[x1,…,xd], which we can write as P(y | x1,…,xd)P(y | x1,…,xd). Bayes's theorem tells us how to express this in terms of quantities we can compute more directly:

P(y | x1,…,xd)=P(x1,…,xd | y)P(y)P(x1,…,xd)P(y | x1,…,xd)=P(x1,…,xd | y)P(y)P(x1,…,xd)

The main assumption in Naive Bayes is that, given the label, the observed features are conditionally independent i.e.

P(x1,…,xd | y)=P(x1 | y)×…×P(xd | y)P(x1,…,xd | y)=P(x1 | y)×…×P(xd | y)

Therefore, we can rewrite Bayes rule as

P(y | x1,…,xd)=P(x1 | y)×…×P(xd | y)P(y)P(x1,…,xd)P(y | x1,…,xd)=P(x1 | y)×…×P(xd | y)P(y)P(x1,…,xd)

Training Naive Bayes¶

One way to train a Naive Bayes classifier is done using frequentist approach to calculate probability, which is simply going over the training data and calculating the frequency of different observations in the training set given different labels. For example,
P(x1=i | y=j)=P(x1=i,y=j)P(y=j)=Number of times in training data x1=i and y=jTotal number of times in training data y=jP(x1=i | y=j)=P(x1=i,y=j)P(y=j)=Number of times in training data x1=i and y=jTotal number of times in training data y=j

Testing Naive Bayes¶

During the testing phase, we try to estimate the probability of a label given an observed feature vector. We combine the probabilities computed from training data to estimate the probability of a given label. For example, if we are trying to decide between two labels y1y1 and y2y2, then we compute the ratio of the posterior probabilities for each label:

P(y1 | x1,…,xd)P(y2 | x1,..,xd)=P(x1,…,xd | y1)P(x1,…,xd | y2)P(y1)P(y2)=P(x1 | y1)×…×P(xd | y1)P(y1)P(x1 | y2)×…×P(xd | y2)P(y2)P(y1 | x1,…,xd)P(y2 | x1,..,xd)=P(x1,…,xd | y1)P(x1,…,xd | y2)P(y1)P(y2)=P(x1 | y1)×…×P(xd | y1)P(y1)P(x1 | y2)×…×P(xd | y2)P(y2)

All we need now is to compute P(x1|yi),…,P(xd | yi)P(x1|yi),…,P(xd | yi) and P(yi)P(yi) for each label by plugging in the numbers we got during training. The label with the higher posterior probabilities is the one that is selected.

4.1 Llama Breed Problem using Naive Bayes [5pts] [W]¶

Above are images of two different breeds of llamas – the Suri and the Wooly. The difference between these two breeds is subtle, as these two breeds are often mixed up. However the Suri Llama is vastly more valuable than the Wooly llama. You devise a way to determine with some confidence, which is which – without the need for expensive genetic testing.

You look at four key features of the llama: {curly hair, over 14 inch tail, over 400 pounds, extremely shy}.

You only have 6 randomly chosen llamas to work with, and their breed as the ground truth. You record the data as vectors with the entry 1 if true and 0 if false. For example a llama with vector {1,1,0,1} would have curly hair, a tail over 14 inches, be less than 400 pounds, and be extremely shy.

The Suri Llamas yield the following data: {1, 0, 1, 0}, {1, 1, 0, 1}, {1, 1, 1, 1}, {0, 0, 0, 1}

The Wooly Llamas yield the following data: {0, 0, 1, 0}, {1, 1, 0, 0}.

Now is the time to test your method!

You see a new llama you are interested in that has curly hair, does not have a tail over 14 inches, is less than 400 pounds, and is not shy.

Using Naive Bayes, is this a Suri or a Wooly Llama?

NOTE: We expect students to show their work (Naive Bayes calculations) and not just the final answer.

Answer: ...

Wooly:

Curly Hair: 1/2 Tail over 14 inches: 1/2 Over 400 pounds = 1/2 Extremly Shy: 0

Suri Llama:

Curly Hair: 3/4 Tail Over 14 Inches: 1/2 Over 400 Pounds: 1/2 Extremly Shy: 1/2

P(Wooly) = 1,0,0,0 = 1/2 4/2 1/2 1 1/3 = 1 / 24 P(Suri) = 1,0,0,0 = 3/4 1/2 1/2 1/2 2/3 = 1 / 16

P(Suri) / P(Wooly) = (1/16) / (1/24) = 1.5

This proves that it is more likely to be a Suri Llama

4.2 News Data Sentiment Classification via Logistic Regression [30pts] **[P]**¶

This dataset contains the sentiments for financial news headlines from the perspective of a retail investor. The sentiment of news has 3 classes, negative, positive and neutral. In this problem, we only use the negative (class label = 0) and positive (class label = 1) classes for binary logistic regression. For data preprocessing, we remove the duplicate headlines and remove the neutral class to get 1967 unique news headlines. Then we randomly split the 1967 headlines into training set and evaluation set with 8:2 ratio. We use the training set to fit a binary logistic regression model.

The code which is provided loads the documents, preprocess the data, builds a “bag of words” representation of each document. Your task is to complete the missing portions of the code in logisticRegression.py to determine whether a news headline is negative or positive.

In logistic_regression.py file, complete the following functions:

  • sigmoid: transform s=xθs=xθ to probability of being positive using sigmoid function, which is 11+e−s11+e−s.
  • bias_augment: augment xx with 1's to account for bias term in θθ
  • predict_probs: predicts the probability of positive label P(y=1|x)P(y=1|x)
  • predict_labels: predicts labels
  • loss: calculates binary cross-entropy loss
  • gradient: calculate the gradient of the loss function with respect to the parameters θθ.
  • accuracy: calculate the accuracy of predictions
  • evaluate: gives loss and accuracy for a given set of points
  • fit: fit the logistic regression model on the training data.
Logistic Regression Overview:¶
  1. In logistic regression, we model the conditional probability using parameters θθ, which includes a bias term b. p(yi=1|xi;θ)=hθ(xi)=σ(xθ)p(yi=1|xi;θ)=hθ(xi)=σ(xθ) p(yi=0|xi;θ)=1−hθ(xi)p(yi=0|xi;θ)=1−hθ(xi)

where σ(⋅)σ(⋅) is the sigmoid function as follows: σ(s)=11+e−sσ(s)=11+e−s

  1. The conditional probabilities of the positive class (y=1)(y=1) and the negative class (y=0)(y=0) of the sample xixi attributes are combined into one equation as follows:
p(yi|xi;θ)=(hθ(xi))yi(1−hθ(xi))1−yip(yi|xi;θ)=(hθ(xi))yi(1−hθ(xi))1−yi
  1. Assuming that the samples are independent of each other, the likelihood of the entire dataset is the product of the probabilities of all samples. We use maximum likehood estimation to estimate the model parameters θθ. The negative log likelihood (scaled by the dataset size NN) is given by: L(θ∣X,Y)=−1NN∑i=1yiloghθ(xi)+(1−yi)log(1−hθ(xi))L(θ∣X,Y)=−1N∑i=1Nyilog⁡hθ(xi)+(1−yi)log⁡(1−hθ(xi))

where:

N=N= number of training samples
xi=xi= bag of words features of the i-th training sample
yi=yi= label of the i-th training sample

Note that this will be our model's loss function

  1. Then calculate the gradient ▽θL▽θL and use gradient descent to optimize the loss function: θt+1=θt−η⋅▽θL(θt∣X,Y)θt+1=θt−η⋅▽θL(θt∣X,Y)

where ηη is the learning rate and the gradient ▽θL▽θL is given by:

▽θL(θ∣X,Y)=1NN∑i=1x⊤i(hθ(xi)−yi)▽θL(θ∣X,Y)=1N∑i=1Nxi⊤(hθ(xi)−yi)

4.2.1 Local Tests for Logistic Regression [No Points]¶

You may test your implementation of the functions contained in logistic_regression.py in the cell below. Feel free to comment out tests for functions that have not been completed yet. See Using the Local Tests for more details.

In [50]:
from utilities.localtests import TestLogisticRegression

unittest_lr = TestLogisticRegression()
unittest_lr.test_sigmoid()
unittest_lr.test_bias_augment()
unittest_lr.test_loss()
unittest_lr.test_predict_probs()
unittest_lr.test_predict_labels()
unittest_lr.test_loss()
unittest_lr.test_accuracy()
unittest_lr.test_evaluate()
unittest_lr.test_fit()
UnitTest passed successfully for "Logistic Regression sigmoid"!
UnitTest passed successfully for "Logistic Regression bias_augment"!
UnitTest passed successfully for "Logistic Regression loss"!
UnitTest passed successfully for "Logistic Regression predict_probs"!
UnitTest passed successfully for "Logistic Regression predict_labels"!
UnitTest passed successfully for "Logistic Regression loss"!
UnitTest passed successfully for "Logistic Regression accuracy"!
UnitTest passed successfully for "Logistic Regression evaluate"!
Epoch 0:
	train loss: 0.675	train acc: 0.7
	val loss:   0.675	val acc:   0.7
UnitTest passed successfully for "Logistic Regression fit"!

4.2.2 Logistic Regression Model Training [No Points]¶

In [51]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from logistic_regression import LogisticRegression as LogReg
In [52]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

news_data = pd.read_csv("./data/news-data.csv",
                    encoding='cp437', header=None)

class_to_label_mappings = {
    "negative": 0,
    "positive": 1
}

label_to_class_mappings = {
    0 : "negative",
    1 : "positive"
}

news_data.columns = ["Sentiment", "News"]
news_data.drop_duplicates(inplace=True)

news_data = news_data[news_data.Sentiment != "neutral"]

news_data["Sentiment"] = news_data["Sentiment"].map(
    class_to_label_mappings)

vectorizer = text.CountVectorizer(stop_words='english')

X = news_data['News'].values
y = news_data['Sentiment'].values.reshape(-1, 1)

RANDOM_SEED = 5
BOW = vectorizer.fit_transform(X).toarray()
indices = np.arange(len(news_data))
X_train, X_test, y_train, y_test, indices_train, indices_test = train_test_split(
    BOW, y, indices, test_size=0.2, random_state=RANDOM_SEED)

Fit the model to the training data. Feel free to try different learning rates lr and number of epochs to achieve >80% test accuracy.

In [53]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

model = LogReg()
lr = 0.05
epochs = 10000
theta = model.fit(X_train, y_train, X_test, y_test, lr, epochs)
Epoch 0:
	train loss: 0.69	train acc: 0.7
	val loss:   0.691	val acc:   0.665
Epoch 1000:
	train loss: 0.436	train acc: 0.794
	val loss:   0.532	val acc:   0.701
Epoch 2000:
	train loss: 0.364	train acc: 0.846
	val loss:   0.484	val acc:   0.746
Epoch 3000:
	train loss: 0.318	train acc: 0.873
	val loss:   0.456	val acc:   0.761
Epoch 4000:
	train loss: 0.286	train acc: 0.896
	val loss:   0.438	val acc:   0.772
Epoch 5000:
	train loss: 0.262	train acc: 0.914
	val loss:   0.425	val acc:   0.782
Epoch 6000:
	train loss: 0.242	train acc: 0.926
	val loss:   0.416	val acc:   0.789
Epoch 7000:
	train loss: 0.226	train acc: 0.933
	val loss:   0.409	val acc:   0.797
Epoch 8000:
	train loss: 0.212	train acc: 0.943
	val loss:   0.404	val acc:   0.802
Epoch 9000:
	train loss: 0.2	train acc: 0.95
	val loss:   0.4	val acc:   0.799

4.2.3 Logistic Regression Model Evaluation [No Points]¶

Evaluate the model on the test dataset

In [54]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

test_loss, test_acc = model.evaluate(X_test, y_test, theta)
print(f"Test Dataset Accuracy: {round(test_acc, 3)}")
Test Dataset Accuracy: 0.807

Plotting the loss function on the training data and the test data for every 100th epoch

In [55]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

model.plot_loss()

Plotting the accuracy function on the training data and the test data for each epoch

In [56]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

model.plot_accuracy()
In [57]:
np.reshape(X_test[0], (1, X_test.shape[1])).shape
Out[57]:
(1, 5286)

Check out sample evaluations from the test set.

In [58]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

num_samples = 10
for i in range(10):
    rand_index = np.random.randint(0, len(X_test))
    x_test = np.reshape(X_test[rand_index], (1, X_test.shape[1]))
    prob = model.predict_probs(model.bias_augment(x_test), theta)
    pred = model.predict_labels(prob)
    print(f"Input News: {X[indices_test[rand_index]]}\n")
    print(f"Predicted Sentiment: {label_to_class_mappings[pred[0][0]]}")
    print(f"Actual Sentiment: {label_to_class_mappings[y_test[rand_index][0]]}\n")
Input News: Earnings per share for January-June 2010 were EUR0 .30 , an increase of 20 % year-on-year EUR0 .25 .

Predicted Sentiment: positive
Actual Sentiment: positive

Input News: Earnings per share ( EPS ) amounted to a loss of EUR0 .38 .

Predicted Sentiment: positive
Actual Sentiment: negative

Input News: Food sales totalled EUR 323.5 mn in October 2009 , representing a decrease of 5.5 % from October 2008 .

Predicted Sentiment: negative
Actual Sentiment: negative

Input News: Kesko pursues a strategy of healthy , focused growth concentrating on sales and services to consumer-customers .

Predicted Sentiment: positive
Actual Sentiment: positive

Input News: Finnair 's Tallinn-based subsidiary , Aero AS , carried 23,335 passengers , a decline of 60.8 % , on routes between Helsinki and the Baltic capitals and within Southern Finland .

Predicted Sentiment: positive
Actual Sentiment: negative

Input News: Finnish pharmaceuticals company Orion 's net sales rose to EUR 190mn in the first quarter of 2009 from EUR 180mn in the first quarter of 2008 .

Predicted Sentiment: positive
Actual Sentiment: positive

Input News: Fiskars , a 360-year-old global business best known for its orange-handled scissors , expects to derive synergies of EUR5 million a year by consolidating certain parts of the housewares division where the two compete .

Predicted Sentiment: positive
Actual Sentiment: positive

Input News: - Net sales for the period are expected to fall well below that of last year and the result after non-recurring items is expected to be in the red .

Predicted Sentiment: negative
Actual Sentiment: negative

Input News: A new production line is being completed for the contract production of hormone treatments .

Predicted Sentiment: positive
Actual Sentiment: positive

Input News: Fiskars , a 360-year-old global business best known for its orange-handled scissors , expects to derive synergies of EUR5 million a year by consolidating certain parts of the housewares division where the two compete .

Predicted Sentiment: positive
Actual Sentiment: positive

Q5: Noise in PCA and Linear Regression [15pts] **[W]**¶

Both PCA and least squares regression can be viewed as algorithms for inferring (linear) relationships among data variables. In this part of the assignment, you will develop some intuition for the differences between these two approaches and develop an understanding of the settings that are better suited to using PCA or better suited to using the least squares fit.

The high level bit is that PCA is useful when there is a set of latent (hidden/underlying) variables, and all the coordinates of your data are linear combinations (plus noise) of those variables. The least squares fit is useful when you have direct access to the independent variables, so any noisy coordinates are linear combinations (plus noise) of known variables.

5.1 Slope Functions [5 pts] **[W]**¶

In the following cell, complete the following:

  1. pca_slope: For this function, assume that XX is the first feature and yy is the second feature for the data. Write a function, that takes in the first feature vector XX and the second feature vector yy. Stack these two feature vectors into a single N x 2 matrix and use this to determine the first principal component vector of this dataset. Be careful of how you are stacking the two vectors. You can check the output by printing it which should help you debug. Finally, return the slope of this first component. You should use the PCA implementation from Q2.

  2. lr_slope: Write a function that takes XX and yy and returns the slope of the least squares fit. You should use the Linear Regression implementation from Q3 but do not use any kind of regularization. Think about how weight could relate to slope.

In later subparts, we consider the case where our data consists of noisy measurements of xx and yy. For each part, we will evaluate the quality of the relationship recovered by PCA, and that recovered by standard least squares regression.

As a reminder, least squares regression minimizes the squared error of the dependent variable from its prediction. Namely, given (xi,yi)(xi,yi) pairs, least squares returns the line l(x)l(x) that minimizes ∑i(yi−l(xi))2∑i(yi−l(xi))2.

In [59]:
import numpy as np
from pca import PCA
from regression import Regression

def pca_slope(X, y):
    """
    Calculates the slope of the first principal component given by PCA

    Args: 
        x: N x 1 array of feature x
        y: N x 1 array of feature y
    Return:
        slope: (float) scalar slope of the first principal component
    """

    raise NotImplementedError



def lr_slope(X, y):
    """
    Calculates the slope of the best fit returned by linear_fit_closed()

    For this function don't use any regularization

    Args: 
        X: N x 1 array corresponding to a dataset
        y: N x 1 array of labels y
    Return:
        slope: (float) slope of the best fit
    """

    raise NotImplementedError

We will consider a simple example with two variables, xx and yy, where the true relationship between the variables is y=4xy=4x. Our goal is to recover this relationship—namely, recover the coefficient “4”. We set X=[0,.02,.04,.06,…,1]X=[0,.02,.04,.06,…,1] and y=4xy=4x. Make sure both functions return 4.

In [60]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
x = np.arange(0, 1.02, 0.02).reshape(-1, 1)

y = 4 * np.arange(0, 1.02, 0.02).reshape(-1, 1)

print("Slope of first principal component", pca_slope(x, y))

print("Slope of best linear fit", lr_slope(x, y))

fig = plt.figure()
plt.scatter(x, y)
plt.xlabel("x")
plt.ylabel("y")

if not STUDENT_VERSION:
    fig.text(0.5, 0.5, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT*0.8)
    
plt.show()
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[60], line 8
      4 x = np.arange(0, 1.02, 0.02).reshape(-1, 1)
      6 y = 4 * np.arange(0, 1.02, 0.02).reshape(-1, 1)
----> 8 print("Slope of first principal component", pca_slope(x, y))
     10 print("Slope of best linear fit", lr_slope(x, y))
     12 fig = plt.figure()

Cell In[59], line 16, in pca_slope(X, y)
      5 def pca_slope(X, y):
      6     """
      7     Calculates the slope of the first principal component given by PCA
      8 
   (...)
     13         slope: (float) scalar slope of the first principal component
     14     """
---> 16     raise NotImplementedError

NotImplementedError: 

5.2 Analysis Setup [5 pts] **[W]**¶

Error in y¶

In this subpart, we consider the setting where our data consists of the actual values of xx, and noisy estimates of yy. Run the following cell to see how the data looks when there is error in yy.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
base = np.arange(0.001, 1.001, 0.001).reshape(-1, 1)
c = 0.5
X = base
y = 4 * base + np.random.normal(loc=[0], scale=c, size=base.shape)

fig = plt.figure()
plt.scatter(X, y)
plt.xlabel("x")
plt.ylabel("y")

if not STUDENT_VERSION:
    fig.text(0.5, 0.5, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)

plt.show()

In following cell, you will implement the addNoise function:

  1. Create a vector XX where X=[x1,x2,...,x1000]=[.001,.002,.003,...,1]X=[x1,x2,...,x1000]=[.001,.002,.003,...,1].

  2. For a given noise level cc, set ^yi∼4xi+N(0,c)=4i/1000+N(0,c)y^i∼4xi+N(0,c)=4i/1000+N(0,c), and ^Y=[^y1,^y2,...,^y1000]Y^=[y^1,y^2,...,y^1000]. You can use the np.random.normal function, where scale is equal to noise level, to add noise to your points.

  3. Notice the parameter x_noise in the addNoise function. When this parameter is set to TrueTrue, you will have to add noise to XX. For a given noise level c, let ^xi∼xi+N(0,c)=i/1000+N(0,c)x^i∼xi+N(0,c)=i/1000+N(0,c), and ^X=[^x1,^x2,....^x1000]X^=[x^1,x^2,....x^1000]

  4. Return the pca_slope and lr_slope values of this XX and ^YY^ dataset you have created where ^YY^ has noise (X=XX=X or ^XX^ depending on the problem).

Hint 1: Refer to the above example on how to add noise to XX or YY

Hint 2: Be careful not to add double noise to XX or YY

In [ ]:
def addNoise(c, x_noise = False, seed = 1):
    """
    Creates a dataset with noise and calculates the slope of the dataset
    using the pca_slope and lr_slope functions implemented in this class.

    Args: 
        c: (float) scalar, a given noise level to be used on Y and/or X
        x_noise: (Boolean) When set to False, X should not have noise added
                 When set to True, X should have noise. 
                 Note that the noise added to X should be different from the 
                 noise added to Y. You should NOT use the same noise you add 
                 to Y here.
        seed: (int) Random seed
    Return:
        pca_slope_value: (float) slope value of dataset created using pca_slope
        lr_slope_value: (float) slope value of dataset created using lr_slope

    """
    np.random.seed(seed) #### DO NOT CHANGE THIS ####
    
    ############# START YOUR CODE BELOW #############
    # TODO: Implement this function
    
    
    ############# END YOUR CODE ABOVE #############
    return pca_slope_value, lr_slope_value

A scatter plot with cc on the horizontal axis and the output of pca_slope and lr_slope on the vertical axis has already been implemented for you.

A sample ^YY^ has been taken for each cc in [0,0.05,0.1,…,0.95,1.0][0,0.05,0.1,…,0.95,1.0]. The output of pca_slope is plotted as a red dot, and the output of lr_slope as a blue dot. This has been repeated 30 times, you can see that we end up with a plot of 1260 dots, in 21 columns of 60, half red and half blue. Note that the plot you get might not look exactly like the TA version and that is fine because you might have randomized the noise slightly differently than how we did it.

NOTE: Here, x_noise = False since we only want Y to have noise.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
pca_slope_values = []
linreg_slope_values = []
c_values = []
s_idx = 0

for i in range(30):
    for c in np.arange(0, 1.05, 0.05):
        
        # Calculate pca_slope_value (psv) and lr_slope_value (lsv)
        psv, lsv = addNoise(c, seed = s_idx)
        
        # Append pca and lr slope values to list for plot function
        pca_slope_values.append(psv)
        linreg_slope_values.append(lsv)
        
        # Append c value to list for plot function
        c_values.append(c)
        
        # Increment random seed index
        s_idx += 1

fig = plt.figure()    
plt.scatter(c_values, pca_slope_values, c='r')
plt.scatter(c_values, linreg_slope_values, c='b')
plt.xlabel("c")
plt.ylabel("slope")

if not STUDENT_VERSION:
    fig.text(0.6, 0.4, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.5, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
    
plt.show()

Error in xx and yy¶

We will now examine the case where our data consists of noisy estimates of both xx and yy. Run the following cell to see how the data looks when there is error in both.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
base = np.arange(0.001, 1, 0.001).reshape(-1, 1)
c = 0.5
X = base + np.random.normal(loc=[0], scale=c, size=base.shape)
y = 4 * base + np.random.normal(loc=[0], scale=c, size=base.shape)

fig = plt.figure()
plt.scatter(X, y)
plt.xlabel("x")
plt.ylabel("y")

if not STUDENT_VERSION:
    fig.text(0.5, 0.5, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.8, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
    
plt.show()

In the below cell, we graph the predicted PCA and LR slopes on the vertical axis against the value of c on the horizontal axis. Note that the graph you get might not look exactly like the TA version and that is fine because you might have randomized the noise slightly differently than how we did it.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
pca_slope_values = []
linreg_slope_values = []
c_values = []
s_idx = 0

for i in range(30):
    for c in np.arange(0, 1.05, 0.05):
        
        # Calculate pca_slope_value (psv) and lr_slope_value (lsv), notice x_noise = True
        psv, lsv = addNoise(c, x_noise = True, seed = s_idx)
        
        # Append pca and lr slope values to list for plot function
        pca_slope_values.append(psv)
        linreg_slope_values.append(lsv)
        
        # Append c value to list for plot function
        c_values.append(c)
        
        # Increment random seed index
        s_idx += 1

fig = plt.figure()
plt.scatter(c_values, pca_slope_values, c='r')
plt.scatter(c_values, linreg_slope_values, c='b')
plt.xlabel("c")
plt.ylabel("slope")

if not STUDENT_VERSION:
    fig.text(0.5, 0.5, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.5, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
    
plt.show()

5.3. Analysis [5 pts] **[W]**¶

Based on your observations from previous subsections answer the following questions about the two cases (error in YY and error in both XX and YY) in 2-3 lines.

NOTE:

  1. The closer the value of slope to actual slope ("4" here) the better the algorithm is performing.
  2. You don't need to provide a mathematical proof for this question.
  3. Understanding how PCA and Linear Regression work should help you decipher which case was better for which algorithm. Base your answer on this understanding of how either algorithms works.

QUESTIONS:

  1. Which case does PCA perform worse in? Why does PCA perform worse in this case? (2 Pts)
  2. Why does PCA perform better in the other case? (1 Pt)
  3. Which case does Linear Regression perform well? Why does Linear Regression perform well in this case? (2 Pts)

ANSWERS:

  1. ... Because PCA attempts ot maximize the avariance and data and there is an error in Y, the direction of the max vairance becomes more aligned with the Y-axis. Therefore, PCA performs the worst in the event there is an error in Y.
  2. ... In the case where there is an error in X, PCA performs better because the relationship between x and y is more aligns the max variance.
  3. ... Linear regression works well when there is an error in both XX and YY because it can still identify the underlying trend. PCA can improve the performance of linear regression by reducing the dimensionality of the data and removing noise and redundancy in the independent variables.

Q6 Feature Reduction Implementation [25pts Bonus for All] **[P]** | **[W]**¶

6.1 Implementation [18 Points] **[P]**¶

Feature selection is an integral aspect of machine learning. It is the process of selecting a subset of relevant features that are to be used as the input for the machine learning task. Feature selection may lead to simpler models for easier interpretation, shorter training times, avoidance of the curse of dimensionality, and better generalization by reducing overfitting.

In the feature_reduction.py file, complete the following functions:

  • forward_selection
  • backward_elimination

These functions should each output a list of features.

Forward Selection:¶

In forward selection, we start with a null model, start fitting the model with one individual feature at a time, and select the feature with the minimum p-value. We continue to do this until we have a set of features where one feature's p-value is less than the confidence level.

Steps to implement it:

  1. Choose a significance level (given to you).
  2. Fit all possible simple regression models by considering one feature at a time.
  3. Select the feature with the lowest p-value.
  4. Fit all possible models with one extra feature added to the previously selected feature(s).
  5. Select the feature with the minimum p-value again. if p_value < significance, go to Step 4. Otherwise, terminate.

Backward Elimination:¶

In backward elimination, we start with a full model, and then remove the insignificant feature with the highest p-value (that is greater than the significance level). We continue to do this until we have a final set of significant features.

Steps to implement it:

  1. Choose a significance level (given to you).
  2. Fit a full model including all the features.
  3. Select the feature with the highest p-value. If (p-value > significance level), go to Step 4, otherwise terminate.
  4. Remove the feature under consideration.
  5. Fit a model without this feature. Repeat entire process from Step 3 onwards.

HINT 1: The p-value is known as the observed significance value for a null hypothesis. In our case, the p-value of a feature is associated with the hypothesis H0:βj=0H0:βj=0. If βj=0βj=0, then this feature contributes no predictive power to our model and should be dropped. We reject the null hypothesis if the p-value is smaller than our significance level. More briefly, a p-value is a measure of how much the given feature significantly represents an observed change. A lower p-value represents higher significance. Some more information about p-values can be found here: https://towardsdatascience.com/what-is-a-p-value-b9e6c207247f

HINT 2: For this function, you will have to install statsmodels if not installed already. To do this, run pip install statsmodels in command line/terminal. In the case that you are using an Anaconda environment, run conda install -c conda-forge statsmodels in the command line/terminal. For more information about installation, refer to https://www.statsmodels.org/stable/install.html. The statsmodels library is a Python module that provides classes and functions for the estimation of many different statistical models, as well as for conducting statistical tests, and statistical data exploration. You will have to use this library to choose a regression model to fit your data against. Some more information about this module can be found here: https://www.statsmodels.org/stable/index.html

HINT 3: For step 2 in each of the forward and backward selection functions, you can use the sm.OLS function as your regression model. Also, do not forget to add a bias to your regression model. A function that may help you is the sm.add_constants function.

TIP 4: You should be able to implement these function using only the libraries provided in the cell below.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from feature_reduction import FeatureReduction
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
bc_dataset = load_breast_cancer()
bc = pd.DataFrame(bc_dataset.data, columns = bc_dataset.feature_names)
bc['Diagnosis'] = bc_dataset.target
X = bc.drop('Diagnosis', 1)
y = bc['Diagnosis']
featureselection = FeatureReduction()
#Run the functions to make sure two lists are generated, one for each method
print("Features selected by forward selection:", FeatureReduction.forward_selection(X, y))
print("Features selected by backward elimination:", FeatureReduction.backward_elimination(X, y))

6.2 Feature Selection - Discussion [7pts] **[W]**¶

Question 6.2.1:¶

We have seen two regression methods namely Lasso and Ridge regression earlier in this assignment. Another extremely important and common use-case of these methods is to perform feature selection. Considering there are no restrictions set on the dataset, according to you, which of these two methods is more appropriate for feature selection generally (choose one method)? Why? (3 pts)

Answer: ...

Lasso regression is more appropriate for feature selection than Ridge regression because it can exclude unimportant features by shrinking their coefficients to zero. However, if all features are known to be important, Ridge regression may be more suitable as it does not eliminate any features.

Question 6.2.2:¶

We have seen that we use different subsets of features to get different regression models. These models depend on the relevant features that we have selected. Using forward selection, what fraction of the total possible models can we explore? Assume that the total number of features that we have at our disposal is NN. Remember that in stepwise feature selection (like forward selection and backward elimination), we always include an intercept in our model, so you only need to consider the NN features. (4 pts)

Answer: ...

There are 2^N possible subsets of features and we are adding one of the remaining features so that the total number of models explored in (N(N+1))/2. The fraction we explore is (N(N+1))/(2*2^N)

Q7: Netflix Movie Recommendation Problem Solved using SVD [10pts Bonus for All] **[P]**¶

Let us try to tackle the famous problem of movie recommendation using just our SVD functions that we have implemented. We are given a table of reviews that 600+ users have provided for close to 10,000 different movies. Our challenge is to predict how much a user would rate a movie that they have not seen (or rated) yet. Once we have these ratings, we would then be able to predict which movies to recommend to that user.

Understanding How SVD Helps in Movie Recommendation¶

We are given a dataset of user-movie ratings (RR) that looks like the following:

Ratings in the matrix range from 1-5. In addition, the matrix contains nan wherever there is no rating provided by the user for the corresponding movie. One simple way to utilize this matrix to predict movie ratings for a given user-movie pair would be to fill in each row / column with the average rating for that row / column. For example: For each movie, if any rating is missing, we could just fill in the average value of all available ratings and expect this to be around the actual / expected rating.

While this may sound like a good approximation, it turns out that by just using SVD we can improve the accuracy of the predicted rating.

How does SVD fit into this picture?¶

Recall how we previously used SVD to compress images by throwing out less important information. We could apply the same idea to our above matrix (RR) to generate another matrix (R_R_) which will provide the same information, i.e ratings for any user-movie pairs but by combining only the most important features.

Let's look at this with an example:

Assume that decomposition of matrix RR looks like:

R=UΣVTR=UΣVT

We can re-write this decomposition as follows:

R=U√Σ√ΣVTR=UΣΣVT

If we were to take only the top K singular values from this matrix, we could again write this as:

R_=U√Σk√ΣkVTR_=UΣkΣkVT

Thus we have now effectively separated our ratings matrix RR into two matrices given by: Uk=U[:k]√ΣkUk=U[:k]Σk and Vk=√ΣkVT[:k]Vk=ΣkV[:k]T

There are many ways to visualize the importance of UU and VV matrices but with respect to our context of movie ratings, we can visualize these matrices as follows:

We can imagine each row of UkUk to be holding some information how much each user likes a particular feature (feature1, feature2, feature 3...feature kk). On the contrary, we can imagine each column of VTkVkT to be holding some information about how much each movie relates to the given features (feature 1, feature 2, feature 3 ... feature kk).

Lets denote the row of UkUk by uiui and the column of VTkVkT by mjmj. Then the dot-product: ui⋅mjui⋅mj can provide us with information on how much a user i likes movie j.

What have we achieved by doing this?¶

Starting with a matrix RR containing very few ratings, we have been able to summarize the sparse matrix of ratings into matrices UkUk and VkVk which each contain feature vectors about the Users and the Movies. Since these feature vectors are summarized from only the most important K features (by our SVD), we can predict any User-Movie rating that is closer to the actual value than just taking any average rating of a row / column (recall our brute force solution discussed above).

Now this method in practice is still not close to the state-of-the-art but for a naive and simple method we have used, we can still build some powerful visualizations as we will see in part 3.

We have divided the task into 3 parts:

1) Implement recommender_svd to return matrices UkUk and VkVk

2) Implement predict to predict top 3 movies a given user would watch

3) (Ungraded) Feel free to run the final cell labeled to see some visualizations of the feature vectors you have generated

Hint: Movie IDs are IDs assigned to the movies in the dataset and can be greater than the number of movies. This is why we have given movies_index and users_index as well that map between the movie IDs and the indices in the ratings matrix. Please make sure to use this as well.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from svd_recommender import SVDRecommender
from regression import Regression
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

recommender = SVDRecommender()
recommender.load_movie_data()
regression = Regression()
# Read the data into the respective train and test dataframes
train, test = recommender.load_ratings_datasets()
print("---------------------------------------------")
print("Train Dataset Stats:")
print("Shape of train dataset: {}".format(train.shape))
print("Number of unique users (train): {}".format(train['userId'].unique().shape[0]))
print("Number of unique users (train): {}".format(train['movieId'].unique().shape[0]))
print("Sample of Train Dataset:")
print("------------------------------------------")
print(train.head())
print("------------------------------------------")
print("Test Dataset Stats:")
print("Shape of test dataset: {}".format(test.shape))
print("Number of unique users (test): {}".format(test['userId'].unique().shape[0]))
print("Number of unique users (test): {}".format(test['movieId'].unique().shape[0]))
print("Sample of Test Dataset:")
print("------------------------------------------")
print(test.head())
print("------------------------------------------")

# We will first convert our dataframe into a matrix of Ratings: R
# R[i][j] will indicate rating for movie:(j) provided by user:(i)
# users_index, movies_index will store the mapping between array indices and actual userId / movieId
R, users_index, movies_index = recommender.create_ratings_matrix(train)
print("Shape of Ratings Matrix (R): {}".format(R.shape))

# Replacing `nan` with average rating given for the movie by all users
# Additionally, zero-centering the array to perform SVD
mask = np.isnan(R)
masked_array = np.ma.masked_array(R, mask)
r_means = np.array(np.mean(masked_array, axis=0))
R_filled = masked_array.filled(r_means)
R_filled = R_filled - r_means

7.1.1 Implement the recommender_svd method to use SVD for Recommendation [5pts] **[P]**¶

In svd_recommender.py file, complete the following function:

  • recommender_svd: Use the above equations to output UkUk and VkVk. You can utilize the svd and compress methods from imgcompression.py to retrieve your UU, ΣΣ and VV matrices.

Local Test for recommender_svd Function [No Points]¶

You may test your implementation of the function in the cell below. See Using the Local Tests for more details.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

from utilities.localtests import TestSVDRecommender

unittest_svd_rec = TestSVDRecommender()
unittest_svd_rec.test_recommender_svd()
In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

# Implement the method `recommender_svd` and run it for the following values of features
no_of_features = [2, 3, 8, 15, 18, 25, 30]
test_errors = []

for k in no_of_features:
    U_k, V_k = recommender.recommender_svd(R_filled, k)
    pred = [] # to store the predicted ratings
    for _, row in test.iterrows():
        user = row['userId']
        movie = row['movieId']
        u_index = users_index[user]
        # If we have a prediction for this movie, use that
        if movie in movies_index:
            m_index = movies_index[movie]
            pred_rating = np.dot(U_k[u_index, :], V_k[:,m_index]) + r_means[m_index]
        # Else, use an average of the users ratings
        else:
            pred_rating = np.mean(np.dot(U_k[u_index], V_k)) + r_means[m_index]
        pred.append(pred_rating)
    test_error = regression.rmse(test['rating'], pred)
    test_errors.append(test_error)
    print("RMSE for k = {} --> {}".format(k, test_error))

Plot the Test Error over the different values of k

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################
fig = plt.figure()
plt.plot(no_of_features, test_errors, 'bo')
plt.plot(no_of_features, test_errors)
plt.xlabel("Value for k")
plt.ylabel("RMSE on Test Dataset")
plt.title("SVD Recommendation Test Error with Different k values")

if not STUDENT_VERSION:
    fig.text(0.5, 0.5, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.5, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
    
plt.show()

7.1.2 Implement the predict method to find which movie a user is interested in watching next [5pts] **[P]**¶

Our goal here is to predict movies that a user would be interested in watching next. Since our dataset contains a large list of movies and our model is very naive, filtering among this huge set for top 3 movies can produce results that we may not correlate immediately. Therefore, we'll restrict this prediction to only movies among a subset as given by movies_pool.

Let us consider a user (ID: 660) who has already watched and rated well (>3) on the following movies:

  • Iron Man (2008)
  • Thor: The Dark World (2013)
  • Avengers, The (2012)

The following cell tries to predict which among the movies given by the list below, the user would be most interested in watching next:
movies_pool:

  • Ant-Man (2015)
  • Iron Man 2 (2010)
  • Avengers: Age of Ultron (2015)
  • Thor (2011)
  • Captain America: The First Avenger (2011)
  • Man of Steel (2013)
  • Star Wars: Episode IV - A New Hope (1977)
  • Ladybird Ladybird (1994)
  • Man of the House (1995)
  • Jungle Book, The (1994)

In svd_recommender.py file, complete the following function:

  • predict: Predict the next 3 movies that the user would be most interested in watching among the ones above.

HINT: You can use the method get_movie_id_by_name to convert movie names into movie IDs and vice-versa.

NOTE: The user may have already watched and rated some of the movies in movies_pool. Remember to filter these out before returning the output. The original Ratings Matrix, RR might come in handy here along with np.isnan

Local Test for predict Functions [No Points]¶

You may test your implementation of the function in the cell below. See Using the Local Tests for more details.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

unittest_svd_rec.test_predict()

7.2 Visualize Movie Vectors [No Points]¶

Our model is still a very naive model, but it can still be used for some powerful analysis such as clustering similar movies together based on user's ratings.

We have said that our matrix VkVk that we have generated above contains information about movies. That is, each column in VkVk contains (feature 1, feature 2, .... feature kk) for each movie. We can also say this in other terms that VkVk gives us a feature vector (of length k) for each movie that we can visualize in a kk-dimensional space. For example, using this feature vector, we can find out which movies are similar or vary.

While we would love to visualize a kk-dimensional space, the constraints of our 2D screen wouldn't really allow us to do so. Instead let us set K=2K=2 and try to plot the feature vectors for just a couple of these movies.

As a fun activity run the following cell to visualize how our model separates the two sets of movies given below.

NOTE: There are 2 possible visualizations. Your plot could be the one that's given on the expected PDF or the one where the y-coordinates are inverted.

In [ ]:
###############################
### DO NOT CHANGE THIS CELL ###
###############################

marvel_movies = ['Thor: The Dark World (2013)',
                'Avengers: Age of Ultron (2015)',
                'Ant-Man (2015)',
                'Iron Man 2 (2010)',
                'Avengers, The (2012)',
                'Thor (2011)',
                'Captain America: The First Avenger (2011)']
marvel_labels = ['Blue'] * len(marvel_movies)
star_wars_movies = [
                'Star Wars: Episode IV - A New Hope (1977)',
                'Star Wars: Episode V - The Empire Strikes Back (1980)',
                'Star Wars: Episode VI - Return of the Jedi (1983)',
                'Star Wars: Episode I - The Phantom Menace (1999)',
                'Star Wars: Episode II - Attack of the Clones (2002)',
                'Star Wars: Episode III - Revenge of the Sith (2005)',
]
star_wars_labels = ['Green'] * len(star_wars_movies)


movie_titles = star_wars_movies + marvel_movies
genre_labels = star_wars_labels + marvel_labels

movie_indices = [movies_index[recommender.get_movie_id_by_name(str(x))] for x in movie_titles]

_, V_k = recommender.recommender_svd(R_filled, k=2)
x, y = V_k[0, movie_indices], V_k[1, movie_indices]
fig = plt.figure()
plt.scatter(x, y, c=genre_labels)
for i, movie_name in enumerate(movie_titles):
    plt.annotate(movie_name, (x[i], y[i]))
    
if not STUDENT_VERSION:
    fig.text(0.5, 0.5, EO_TEXT, transform=fig.transFigure,
        fontsize=EO_SIZE/2, color=EO_COLOR, alpha=EO_ALPHA*0.5, fontname=EO_FONT,
        ha='center', va='center', rotation=EO_ROT)
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]: